AFRL-SR-BL-TR-99- 


1999 


report  documentation  pa 


WOM  f»00rtMiq  Wfowi  '!UiSSSlImS«  ^MW^iufTin  tenicn.  Oirwnrm  for  imonutioon  OMnoom  aM  M 

"l.  AfifNCV  USI  OHU  {t*«*  oiMUtt  1 2*  "‘i^May  99  I _ Final  Technical  I 

TrrTL^MirsuiTmJ  I 

«•  IIIW  MM  «««  _  -rxrxjr-r  a  m  m  T?  I  c 


mn^  MU  lourcfft. 
otMT  noact  ot  wi» 
nm.  1219  Jtf^Miofi 
2QM3. 


TuU  AND 

VISCOUS  EFFECTS  IN  PORE  SCALE  MODELING  OF  IMMISCIBLE 
VISCOUS  _ _ 


r«.  AUTHOmsi 


Durnford,  Deanna  S. 
Trantham,  Heather  Hevelone 


1  WtfOWMtNg  OROAXIZATIOW  WAMt(S)  *NO 
I.PWOnmns,  University 

Fort  Collin3>  CO  80523 


S.  FUNDING  NUMSCKS 


F49620-95-1-0413 


t.  PtMOMMINQ  ORGANIZATION 
UPORT  NUMRCII 


^^^^55s5«Jio7MoS^^ 


Bolling,  AFB 


lo.  SPONSORING  /  MONiTORiNQ 
AGENOr  RERORT  NUMRER 


[t1.  SUPfUMEHTARY  MOTM 


19990803  087 


'tli^RmibmTCMENTA 

Approved  for  Public  Release 
Distribution  Unlimited 


rtEb.  WSTRISUTION  CODE 


13.  AiSTRACT  (MAJomomiOOwofdrt  i!niiiH<!  ( DNAPLs")  such  as  chlorinated  solvents  are, 

Pore  scale  models  capture  the  essential  physics  of  the 

to  solve  the  large  systems  of  equations  inherent  tn  be  usrf  »  predict  source  zone 

labLtory  experiments.  Model  and  laboratory  results  ivere  found  to  compare  well  . 

14.  SURIICT  TERMS  — 

_ _ t- ay  r'cint-flmlnation,  mddeling,  DNAPLs,  MS^RK^COOI^^^^ 


groundwater  contamination,  mddeling,  DNAPLs,  14.  PRICE  CODE 

1  chlorinated  solvents  _  -  . . .  ■- 

I  _ _ LHiiULIIOir^  20.  UMITATIOH  OP  AISTRACT 

“c’l^ifled  .  bucussifled _ \ _ Unclassified 


NSN  7540-01 -2a0-5500 


Standard  Forni  298  {Rdv.  2-89) 

Pt giCHBM  AMSI  S«.  a**U 


EXECUTIVE  SUMMARY 


Deanna  Dumford 
Department  of  Civil  Engineering 
Colorado  State  University 

Heather  Hevelone  Trantham 
Department  of  Civil  Engineering 
Colorado  State  University 

Dense  Non- Aqueous  Phase  Liquid  (DNAPL)  source  zone  characterization  is 
important  for  risk  assessments,  feasibility  studies  and  identification  of  appropriate 
remediation  technologies  at  DNAPL-contaminated  sites.  The  extent  and  configuration  of 
the  source  zone  are  also  important  input  to  multiphase  mass  transfer  models.  When  a 
DNAPL  displaces  water,  the  interface  between  fluids  is  often  unstable  due  to  fluid  and 
porous  media  properties.  Solutions  stemming  from  the  continuity  equation  do  not 
consider  the  stability  of  the  interface  at  the  pore  scale.  However,  many  relevant  questions 
need  to  be  answered  by  pore  scale  models.  How  much  lateral  spreading  of  the  DNAPL 
can  be  expected?  How  much  of  the  pore  space  was  affected  by  the  DNAPL?  How  much 
DNAPL  is  in  contact  with  the  groundwater?  The  model  discussed  in  this  abstract  can 
answer  these  questions. 

The  stochastic  aggregation  model  (Trantham  and  Durnford,  1999)  uses  a  modified 
DLA  algorithm  (Witten  and  Sander,  1983)  to  model  the  displacement  of  water  by 
DNAPLs.  Using  the  model,  it  is  possible  to  look  at  the  displacement  of  these  fluids  under 
different  flow  rates,  changing  porous  media  properties,  and  in  a  gravity  field.  The  fi’ont 
configurations  produced  under  these  conditions  range  from  viscous  fingers  to  flat  fi'onts. 
The  model  uses  macroscopic,  dimensionless  capillary  and  bond  numbers  to  circumvent 
the  difficulty  of  defining  parameters  on  the  microscale. 

A  version  of  the  model  is  available  for  use  on  the  www.colostate.edu/Depts/CE/ 
(click  on  “what’s  happening”).  Trantham  and  Dumford  (1999)  can  also  be  downloaded 
from  this  site. 

Equation  [1]  is  the  non-dimensionalized  capillary  pressure  gradient.  This  is 
equivalent  to  the  stability  criteria  used  by  Saffinan  and  Taylor  (1958)  and  Chuoke  (1959). 
The  equation  has  been  macroscopically  non-dimensionalized  following  Hilfer  and  0ren 
(1996): 

Ap  _  A/g/sing 

Ic  P  1r  P  P 


The  macroscopic  flux  is  given  by  q,  (x  is  the  viscosity,  k  is  the  intrinsic 
permeability.  Pa  is  the  displacement  pressure,  1  is  a  macroscopic  length,  and  a  is  a  scaling 
factor.  Subscripts  w  and  nw  refer  to  the  wetting  and  non-wetting  phases,  respectively. 
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The  scaling  factor  is  necessary  when  using  a  macroscopic  equation  to  describe  the  pore- 
scale  movement  of  the  interface.  The  density  difference  between  the  non-wetting  and 
wetting  phases  is  given  by  Ap,  and  sin  a  accounts  for  any  angle  of  inclination  from  the 
horizontal. 

The  first  two  terms  in  Equation  [1]  are  the  macroscopic  capillary  numbers  for  the 
wetting  and  non-wetting  phases,  respectively.  The  last  term  in  Equation  [1]  is  the 
macroscopic  bond  number.  The  linear  combination  of  bond  and  capillary  numbers  is 
called  the  transition  number.  The  transition  number  quantifies  the  stability  of  the 
DNAPL-water  front  as  it  transitions  from  viscous  fingers  to  a  flat  front.  The  transition 
number  is  given  as 


r  = 


[2] 


In  the  stochastic  aggregation  model,  the  sticking  probability  (Ps),  accounts  for  the 
changes  in  front  stability.  A  unique  relationship  between  the  sticking  probability  and  the 
transition  number  was  developed  from  43  two-dimensional,  DNAPL-water  laboratory 
experiments.  The  relationship  between  the  sticking  probability  and  the  transition  number 
was  determined  as 


P^  =  0.1227822|rp'-^^^^'’  [3] 


The  model  is  validated  for  homogeneous  porous  media  through  comparison  with 
laboratory  experiments  using  Tetrachloroethylene  (PERC),  1,2  Dichloroethane  (DCA), 
Carbon  Tetrachloride  (CTET)  and  Mobile  Pyrogard  53.  Examples  are  shown  in  Figure  1. 
The  laboratory  experiments  are  shown  on  the  left  in  Figure  1  with  the  corresponding 
model  simulation  shown  on  the  right.  It  is  important  to  compare  the  characteristics  of  the 
DNAPL-water  displacement,  i.e.  branching,  spreading,  and  relative  thickness  of  the 
fingers. 
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Figure  1.  Application  of  model  to  DNAPL-water  displacement  in  a  two-dimensional  flume,  (a)  PERC  90°, 
0=3  ml  min*',  min,  T=0.354,  Ps=0.485.  (b)  1,2  Dichloroethane  45°,  Q=6  ml  min’',  t=l  min,  T=0.473, 
Ps=0.331.  (c)  Carbon  Tetrachloride  30°,  Q=6  ml  min*',  t=l  min,  T=0.194,  Ps=1.0.  (d)  Mobile  Pyrogard  53, 

90°,  0=3  ml  min’’,  t=l  min,  T=-2.44,  Ps=0.037. 

The  model  is  validated  for  heterogeneous  porous  media  through  comparison  with 
the  lab  experiment  of  Kueper  and  Frind  (1991).  Figure  2  shows  the  heterogeneous  sand 
packing  used  by  Kueper  and  Frind  (1991).  PERC  was  allowed  to  infiltrate  from  the 
source  area  at  the  top  of  the  flume  under  a  constant  head  of  4  cm.  Figure  3  shows  the 
results  of  Kueper  and  Frind’s  (1991)  lab  experiment  on  the  right  with  the  corresponding 
model  simulation  on  the  left.  The  model  simulations  capture  the  characteristic  DNAPL- 
water  displacement  and  show  that  the  PERC  preferentially  travels  through  regions  of 
higher  permeability. 
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Figure  2.  The  heterogeneous  sand  packing  of  Kueper  and  Frind’s  (1991)  PERC-water 
displacement  experiment. 
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Figure  3.  Model  simulations  (on  the  left)  compared  to  the  experimental  results  of  Kueper  and  Frind  (1991) 
shown  on  the  right,  (a)  34  seconds,  (b)  126  seconds,  (c)  245  seconds. 


The  model  is  applied  to  the  determination  of  the  bulk  retention  capacity.  Bulk 
retention  capacity  is  a  measure  of  the  amount  of  DNAPL  retained  per  unit  volume  of 
aquifer  (Johnson  and  Kueper,  1996).  It  is  used  in  the  determination  of  DNAPL  depth. 
Changes  in  bulk  retention  capacity  are  investigated  with  changing  inflow  rate,  intrinsic 
permeability,  and  DNAPL  viscosity. 

The  model  is  currently  well  formulated  for  simulating  two-dimensional  DNAPL- 
water  displacement  for  a  variety  of  DNAPLs  under  the  conditions  of  different  inflow 
rates,  different  angles  of  inclination  from  the  horizontal,  and  varying  heterogeneous 
media  properties.  The  comparison  of  model  simulations  to  laboratory  experiments 
remains  limited  to  visual  inspection,  as  is  the  case  for  other  models  currently  attempting 
to  model  DNAPL-water  displacements.  We  are  currently  working  on  expanding  the 
model  to  three  dimensions  and  exploring  better  means  of  graphical  representation  of  the 
model’s  output. 

The  dissertation  submitted  by  Heather  Trantham-Hevelone  to  Colorado  State 
University  is  attached  as  an  appendix  to  this  report,  without  appendices.  These  are 
available  from  the  author  by  request. 
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CHAPTER  1 


INTRODUCTION 

Groundwater  contamination  by  non-aqueous  phase  liquids  (NAPLs)  has  received 
increasing  attention  in  the  last  decade.  Light  Non-Aqueous  Phase  Liquids  (LNAPLs)  are 
liquids  with  specific  gravities  less  than  water.  They  are  associated  with  the  production, 
refining,  and  distribution  of  petroleum  products  (Bedient,  1994).  Dense  Non- Aqueous 
Phase  Liquids  (DNAPLs)  are  liquids  with  specific  gravities  greater  than  water.  Many 
commonly  used  DNAPLs  are  chlorinated  solvents,  but  there  are  several  that  are 
halogenated  organics,  PCB  mixtures,  and  pesticides.  Typical  uses  for  DNAPLs  include 
metal  cleaning  and  degreasing,  dry  cleaning,  paint  removal,  and  use  as  adhesives  and 
aerosols  (Pankow  et  al.,  1996).  DNAPL  solubilities  in  air  and  water  are  high  enough  to 
far  exceed  regulatory  standards  when  dissolution  occurs  in  the  groundwater  even  though 
they  are  essentially  immiscible  in  water  (Mabey  et  al.,  1982).  Both  lighter  and  denser 
NAPLs  can  contaminate  water  supply  systems,  but  due  to  unfavorable  viscosity  and 
density  differences  with  water,  DNAPLs  can  penetrate  the  water  table  and  migrate 
downward.  As  vertical  migration  occurs,  a  trail  of  trapped  residual  is  left  behind.  If 
lower  permeability  soils  are  encountered,  small  pools  of  DNAPL  may  form.  Both 
scenarios  can  lead  to  long  term  contamination. 

DNAPL  source  zone  characterization  is  important  for  risk  assessments,  feasibility 
studies  and  identification  of  appropriate  remediation  technologies  at  DNAPL- 


contaminated  sites.  The  extent  and  configuration  of  the  source  zone  are  also  important 
input  to  multiphase  mass  transfer  models.  Field  and  laboratory  experiments  to  determine 
source  zones  are  either  not  feasible  or  not  permitted  due  to  the  hazardous  nature  of 
DNAPLs.  Numerical  modeling  of  DNAPL  migration  provides  a  solution  to  this  problem, 
however  accurate  qualitative  and  quantitative  descriptions  of  DNAPL  source  zones  are 
unattainable  with  commonly  used  modeling  techniques. 

In  this  study,  the  problem  of  accurate  source  zone  characterization  is  looked  at  in 
terms  of  the  DNAPL-water  interface  stability  on  the  pore  scale.  The  objectives  of  this 
study  are: 

•  to  develop  a  modeling  technique  for  simulating  DNAPL-water  displacements  at 
the  pore  scale. 

•  to  evaluate  the  modeling  technique  through  comparisons  with  DNAPL-water 
displacement  studies  in  homogeneous  porous  media. 

•  to  apply  the  new  modeling  technique  to  DNAPL-water  displacements  in 
heterogeneous  media. 

•  to  investigate  the  use  of  the  modeling  technique  in  calculating  bulk  retention 
capacities. 

Continuum  models  are  used  to  simulate  the  immiscible  transport  of  NAPLs 
through  water-saturated  media.  However,  there  is  no  longer  any  question  that  small  scale 
heterogeneities  control  the  movement  of  DNAPL  through  soil  (Schwille,  1988;  Kueper 
and  Frind,  1991a  and  b;  Kueper  and  Gerhand,  1995,  Illangakesare  et  al.,  1995;  Held  and 
Illangasekare,  1995a  and  b;  and  Illangasekare  et  al.,  1996).  While  these  heterogeneities 
can  be  characterized  statistically,  the  commonly  available  models  are  ill  equipped  to 
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accommodate  their  influence  on  multiphase  transport  (Mercer  et  ah,  1996).  Using 
continuum  models  to  determine  fluid-fluid  displacements  at  the  pore  scale  is 
computationally  prohibitive  and  beyond  the  capability  of  available  numerical  simulators 
(Poulsen  and  Kueper,  1992).  Invasion  percolation  (IP)  models  are  used  for  pore  scale 
fluid-fluid  displacement  (Chatzis  and  Dullien,  1985;  Li  et  al,  1986;  Wardlaw  et  al,  1987; 
Constantinides  and  Payatakes,  1989;  Jerauld  and  Salter,  1990;  Tsakiroglou  and 
Payatakes,  1990;  Ferrand  and  Celia,  1992;  Soil  and  Celia,  1993;  Lowry  and  Miller,  1995; 
Glass  et  al.,  1995).  These  are  applicable  only  for  capillary-dominated,  stable  flow.  IP 
models  neglect  important  processes  because  gravity  and  viscous  forces,  in  addition  to 
capillary  forces  govern  the  flow  patterns  of  many  DNAPLs  of  interest  in  contaminant 
hydrology  (Schwille,  1988). 

DNAPL-water  displacements  are  typically  imstable  depending  on  the  fluid 
properties  of  the  DNAPL  relative  to  water.  Witten  and  Sander’s  (1983)  diffusion  limited 
aggregation  (DLA)  technique  has  been  used  to  model  such  unstable  processes  as 
dielectric  breakdown  (Niemeyer  et  al.,  1984;  Wiesmann  and  Zeller,  1986),  diffusion- 
controlled  polymerization  (Kaufinan  et  al.,  1986),  chemical  dissolution  processes 
(Daccord  et  al.,  1986),  solute  leaching  (Flury  and  Fluhler,  1995),  viscous  fingering 
(Maloy  et  al.,  1985;  Chen  and  Wilkinson,  1985),  fluid-fluid  displacement  in  Hele  Shaw 
cells  (Kadanoff,  1985;  Liang,  1986);  and  fluid-fluid  flow  in  porous  media  in  the  absence 
of  gravity  (Paterson,  1984;  Lenormand,  1987;  Lenormand  et  al.,  1988;  Kiriakidis  et  al., 
1991;  Fernandez  et  al.,  1991). 

The  stochastic  aggregation  modeling  algorithm,  developed  in  this  study,  is  based 
on  Witten  and  Sander’s  (1983)  diffusion  limited  aggregation  (DLA)  technique.  The 
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model  inputs  are  the  essential  properties  governing  front  stability:  DNAPL-water 
viscosity  difference,  DNAPL-water  density  difference,  intrinsic  permeability  of  the 
porous  media,  flow  rate,  and  the  inclination  angle  of  the  porous  media  from  the 
horizontal  if  the  domain  is  two-dimensional.  The  front  stability  is  evaluated  using  criteria 
similar  to  the  instability  criteria  of  SafBnan  and  Taylor  (1958)  and  Chuoke  et  al.  (1959). 
The  background  information  pertaining  to  fluid-fluid  displacements  and  stochastic 
modeling  techniques  is  provided  in  Chapter  2.  The  development  of  the  modeling 
technique  is  presented  in  Chapter  3.  In  Chapter  4,  the  model  is  evaluated  by  comparing 
model  simulations  to  two-dimensional  DNAPL-water  displacements  in  porous  media.  In 
Chapter  5,  the  model  is  used  to  simulate  three  heterogeneous  laboratory  experiments.  An 
application  of  the  model  to  the  evaluation  of  bulk  retention  capacities  is  examined  in 
Chapter  6,  and  the  benefits  and  limitations  of  the  modeling  technique  are  discussed  in 
Chapter  7. 
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CHAPTER  2 


BACKGROUND 

2.1  Fluid-Fluid  Displacement 

Depending  on  the  relative  fluid  properties,  immiscible  fluid  displacements  may 
exhibit  unstable  flow.  An  inherently  unstable  flow  regime  results  when  water  is 
displaced  by  a  DNAPL  with  a  lower  relative  viscosity.  Many  DNAPLs,  including 
chlorinated  solvents,  are  both  more  dense  and  less  viscous  than  water.  The  term 
fingering  is  used  as  a  visual  description  of  these  displacements  in  porous  media  (Chuoke 
et  al.,  1959),  but  does  not  necessarily  imply  unstable  flow.  Unstable  fluid-fluid 
displacements  result  in  fingered  fronts  when  a  small  perturbation  is  caused  by  a 
microscopic  heterogeneity  in  a  macroscopically  homogeneous  media  (Kueper  and  Frind, 
1988).  The  velocity  of  the  displacing  fluid  can  contribute  to  unstable  flow  regimes  (Hill, 
1952)  in  addition  to  the  relative  fluid  properties.  Chuoke  et  al.  (1959)  determined  the 
critical  wavelength  (tip  to  tip  separation  between  the  fingers)  as  a  constraint  accounting 
for  the  stabilizing  or  destabilizing  influence  of  the  surface  tension.  The  problem  is 
compounded  at  the  pore-scale  if  such  details  as  wetting  films,  changing  contact  angles, 
and  heat  and  mass  transfer  are  included  in  the  problem  description. 

The  description  of  fingering  in  homogeneous  porous  media  depends  on  the  length 
scale  chosen,  even  when  small-scale  heterogeneities  and  gravity  are  neglected  (Homsy, 
1987).  Chuoke  et  al.  (1959)  showed  variations  in  fingering  length  scales  with  both 
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increasing  velocity  and  viscosity  difference  in  homogeneous  porous  media.  Applieability 
of  maeroscopic  scales  decreases  with  increasing  velocity,  increasing  mobility,  or 
decreasing  surface  tension  (Homsy,  1987).  Wetting  properties  of  the  two  fluids  are  also 
important.  If  the  invading  fluid  wets  the  medium,  then  the  fingering  can  be  charaeterized 
on  a  macroscopic  continuum  scale.  If  the  invading  fluid  is  non-wetting  with  respect  to 
the  medium,  fingering  is  confined  to  the  pore  scale  (Homsy,  1987). 

There  is  no  preferred  method  of  numerically  modeling  immiscible  instabilities 
(Kueper  and  Frind,  1988).  Analytieal  solutions  require  many  simplifying  assumptions 
due  to  the  nonlinearity  of  the  problem  and  do  not  explicitly  consider  the  stability  of  the 
interface  between  the  fluids.  Other  modeling  techniques  including  stochastic  aggregation 
models  are  equipped  to  model  the  effects  of  pore-scale  instabilities  by  basing  the 
movement  on  decisions  made  at  the  pore  scale. 

2.2  Stochastic  Modeling  Techniques 

The  displacement  of  one  fluid  by  another  in  porous  media  can  be  represented  by 
partiele  aggregation  where  the  growing  aggregate  represents  the  invading  fluid  and  the 
space  surrounding  the  aggregate  represents  the  defending  fluid.  Aggregate  growth  by 
irreversible  particle  collection  is  a  common  phenomenon  in  nature  (Sander,  1984),  and 
several  modeling  techniques  can  be  used  for  partiele  aggregation.  These  include  Invasion 
Percolation  (IP)  modeling  (Wilkinson  and  Willemsen,  1983),  Eden  modeling  (Sander 
1984),  ballistic  aggregation  (Sander,  1984),  and  DLA  (Witten  and  Sander,  1983).  IP 
modeling  is  applicable  to  immiscible  fluid  displacements  dominated  by  capillary  forces 
(Chatzis  and  Dullien,  1985;  Li  et  al.,  1986;  Wardlaw  et  al.,  1987;  Constantinides  and 
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Payatakes,  1989;  Jerauld  and  Salter,  1990;  Tsakiroglou  and  Payatakes,  1990;  Ferrand  and 
Celia,  1992;  Soil  and  Celia,  1993;  Lowry  and  Miller,  1995).  Whereas,  IP  models 
simulate  stable  migration  patterns,  many  DNAPLs  move  in  inherently  imstable  flow 
regimes  (Schwille,  1988).  Other  aggregate  modeling  options  are  needed  to  model  fluid- 
fluid  displacements  where  any  combination  of  viscous,  capillary,  or  gravity  forces  govern 
the  stability  of  the  interface.  Neither  Eden  modeling  nor  ballistic  aggregation  describe 
the  growth  instabilities  characteristic  of  many  processes  including  fluid-fluid 
displacement  with  unfavorable  mobility  ratios  (Sander,  1986).  DLA  models,  however, 
produce  aggregates  resulting  from  unstable  phenomenon  using  a  stochastic  process. 
Following  Witten  and  Sander  (1983),  numerous  other  researchers  have  expanded  on  both 
the  modeling  technique  and  the  applications  for  the  model  (Meakin  and  Deutch,  1986; 
Meakin,  1986;  and  Jullien  et  al.,  1984).  Meakin  (1988)  provides  a  thorough  review  of 
DLA  as  a  model  for  kinetic  growth  phenomena.  A  review  of  IP  and  DLA  modeling  is 
available  in  Appendix  A. 

Lenormand  et  al.  (1988)  performed  physical  micromodel  experiments  and  showed 
dramatically  different  front  configurations  occurring  when  one  fluid  displaces  another  in 
porous  media.  He  divided  the  front  configurations  into  pure  regimes  that  are  modeled  by 
different  stochastic  techniques.  The  domains  of  the  stochastic  modeling  techniques  are 
shown  in  the  phase  diagram  in  Figure  1  (Lenormand,  1987).  The  phase  diagram  is 
applicable  to  displacement  of  the  wetting  fluid  by  the  non-wetting  fluid  without  gravity 
forces.  The  pure  regimes,  labeled  DLA,  IP,  and  Anti-DLA,  correspond  to  displacement 
of  a  viscous  fluid  by  a  non- viscous  fluid,  capillary  displacement,  and  displacement  of  a 
non-viscous  fluid  by  a  viscous  fluid,  respectively.  Anti-DLA  results  in  a  flat  interface 
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and  occurs  when  the  pressure  gradient  is  decreasing  in  the  direction  of  the  flow.  The 
areas  between  these  pure  domains  are  governed  by  a  combination  of  the  applicable 
forces.  The  limits  of  these  domains  are  determined  by  an  evaluation  of  the  capillary 
pressure  at  the  interface  between  the  two  fluids  and  use  the  mobility  ratio  (M)  and  the 
microscopic  capillary  number  (Ca)  as  parameters.  The  mobility  ratio,  M,  is  defined  as 
the  ratio  of  the  dynamic  viscosity  of  the  displacing  fluid  and  the  dynamic  viscosity  of  the 
displaced  fluid.  The  microscopic  capillary  number,  Ca,  is  defined  as  the  ratio  between 
the  viscous  forces  and  the  capillary  forces  (Lake,  1989). 

The  phase  diagram  designates  modeling  techniques  for  the  pure  regimes  of 
capillary  and  viscous  fingering  and  their  limits.  Continuum  models  are  in  the  anti-DLA 
regime.  Most  DNAPLs  of  practical  importance  are  in  the  transition  area  based  on  their 
viscosity,  not  in  one  of  the  pure  domains  shown  in  Figure  2.1.  The  combined  forces 
acting  to  stabilize  or  destabilize  the  front  must  be  considered  to  model  fluid-fluid 
displacements  between  the  regimes.  The  effects  of  gravity  are  also  important  (Wilkinson, 
1984;  Birovljev  et  al.,  1991)  and  should  be  included  in  determining  the  NAPL-water 
front  configuration. 
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Figure  2.1.  The  stochastic  modeling  techniques  corresponding  to  the  regimes 
for  which  their  use  is  suggested  is  shown  with  varying  microscopic  capillary  number  and 
mobility  ratio.  Adapted  from  Lenormand’s  (1987)  phase  diagram  for  immiscible  fluid-fluid  displacement 

in  the  absence  of  gravity. 


Several  studies  have  investigated  the  use  of  DLA  for  modeling  the  displacment  of 
a  more  viscous  fluid  by  a  less  viscous  one  (Paterson,  1984;  Chen  and  Wilkinson,  1985; 
Maloy  et  al.,  1985;  Fernandez  et  al.,  1991;  Kiriakidis,  et  al.,  1991)  in  addition  to 
Lenormand  et  al.  (1988).  Paterson  (1984)  was  the  first  to  suggest  the  application  of  DLA 
to  fluid-fluid  displacements  in  the  absence  of  gravity.  Following  Paterson  (1984),  Maloy 
et  al.  (1985)  conducted  experiments  whereby  air  displaced  epoxy  in  a  radial,  two- 
dimensional  porous  medium.  The  fractal  dimension  of  the  experimental,  tree-like 
structure  was  the  same  as  that  of  DLA  aggregates.  This  study  showed  that  DLA  could  be 


9 


used  to  simulate  the  displacement  of  one  fluid  by  another  in  the  limit  of  an  inviscid 
displacing  fluid.  In  an  effort  to  simulate  a  broader  range  of  displacing  fluids,  Kiriakidis 
et  al.  (1991)  incorporated  a  sticking  probability  into  the  DLA  algorithm.  Kiriakidis  et  al. 
(1991)  based  the  value  of  the  sticking  probability  on  the  pressure  gradient  as  calculated 
by  Lenormand  et  al.  (1988)  and  applied  a  different  sticking  probability  to  each  particle 
depending  on  whether  it  represented  the  wetting  or  non-wetting  fluid.  The  model 
simulations  of  Kiriakidis  et  al.  (1991)  are  not  compared  with  laboratory  experiments. 

One  reason  for  this  may  be  that  the  dynamic  viscosities  corresponding  to  the  simulations 
are  either  much  lower  (7.6  x  10'^  cp)  or  higher  (13  cp)  than  the  majority  of  NAPLs  of 
interest  in  contaminant  hydrology.  With  the  exception  of  PCB’s,  the  majority  of  dynamic 
viscosities  for  NAPLs  lie  in  the  range  of  0.3  cp  to  8  cp  (Mercer  and  Cohen,  1990). 
Fernandez  et  al.  (1991)  also  changed  the  probability  with  which  a  particle  was 
incorporated  into  the  aggregate.  In  this  study,  the  “releasing”  probability  is 
approximately  related  to  the  inverse  of  the  microscopic  capillary  number  through  a 
manipulation  of  the  ratio  of  pressures  (the  pressure  in  the  wetting  or  non-wetting  fluid  as 
it  relates  to  the  capillary  pressure)  and  an  application  of  Darcy’s  law  in  the  absence  of 
gravity.  As  in  the  case  of  Kiriakidis  et  al  (1991),  the  simulations  appear  to  be 
representative  of  NAPL-water  displacements.  It  is  difficult  to  evaluate  the  correctness  of 
this  work  due  to  the  approximate  relationship  between  the  releasing  probability  and  the 
microscopic  capillary  number.  No  evaluation  of  the  microscopic  capillary  number  or  any 
physical  fluid  properties  are  given,  and  hence  it  is  uncertain  if  realistic  fluid  pairs  will  be 
represented  by  this  determination.  In  addition,  Fernandez  et  al.  (1991)  neglects  the 
influence  of  fluid  viscosities  (mobility  ratio). 
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CHAPTERS 


MODEL  FORMULATION  AND  CALIBRATION 

3.1  Model  Formulation 

The  stochastic  aggregation  model  (SAM)  uses  a  modified  DLA  algorithm  to 
model  the  displacement  of  water  by  DNAPLs  which  are  more  or  less  viscous  than  water. 
Using  the  model,  it  is  possible  to  look  at  the  displacement  of  these  fluids  under  conditions 
of  increasing  flow  rates  and  in  a  gravity  field.  The  front  configurations  produced  under 
these  conditions  range  from  viscous  fingers  to  flat  fronts.  The  model  uses  macroscopic, 
dimensionless  capillary  and  bond  numbers  to  circumvent  the  difficulty  of  defining 
parameters  on  the  microscale.  An  analysis  of  the  stability  of  the  front  is  followed  by  an 
explanation  of  the  stochastic  aggregation  algorithm  and  its  relationship  to  the  front 
stability. 

3.1.1  Analysis  of  Front  Stability 

Instabilities  in  the  fluid-fluid  front  are  likely  to  arise  on  the  pore  scale  when  a 
non- wetting  fluid  invades  a  wetting  fluid  in  a  homogeneous  porous  medium  (Homsy, 
1987).  Consider  a  non- wetting  fluid  displacing  a  wetting  fluid  vertically  downwards 
through  a  homogeneous  porous  medium.  This  is  shown  in  Figure  3.1 .  To  initiate  flow 
instability,  the  interface  moves  at  a  pore  scale  velocity,  v,  through  a  distance  dz. 
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Figure  3.1.  Microscopic  interface  perturbation  (adapted  from  Kueper  and  Frind,  1988). 

The  pressures  of  the  wetting  phase  and  non- wetting  phases  at  point  2  are: 

+  [1] 

(.P.),=(P.\*P,gdz-^  [2] 

where  (Pnw)b  (Pnw)2,  (Pw)band  (Pw)2  are  the  pressures  in  the  non- wetting  and  wetting 
phases  at  points  1  and  2.  The  density  and  viscosity  for  each  phase  are  given  by  p  and  p, 
respectively.  The  third  term  on  the  right  hand  sides  of  Equations  [1]  and  [2]  accounts  for 
the  viscous  losses  as  the  interface  moves  from  point  1  to  point  2.  The  k  in  each  equation 
depends  on  the  characteristic  length  dimension  describing  the  cross  section  and  the  shape 
of  the  selected  dimension  representing  the  pore  (Corey,  1994). 

The  change  in  capillary  pressure  between  point  1  and  point  2,  due  to  the 
movement  of  the  interface  is  obtained  by  subtracting  Equation  [2]  from  Equation  [1]  and 
rearranging. 

dPc  =  -  P.)gdz  +  (^  -  ^)vdz  1-3] 

If  the  net  capillary  pressure  is  positive,  then  the  small  displacement  dz  will 
continue  to  propagate.  This  leads  to  a  definition  of  unstable  flow  for  multiphase  flow. 
The  flow  is  unstable  if  the  change  in  capillary  pressure  due  to  displacement  dz  is  positive. 
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This  is  analogous  to  the  analysis  of  Saftman  and  Taylor  (1958)  and  Chuoke  et  al.  (1959) 
where  the  stability  of  the  flow  is  affected  by  both  the  density  difference  and  viscosity 
difference  between  the  two  fluids. 

The  capillary  pressure  gradient  is  given  in  Equation  [4]. 

Vi>.  =  (A.  w 

A  macroscale  equation  similar  to  Equation  [4]  is  derived  starting  with  Darcy’s 
law.  In  immiscible  fluid  flow,  the  capillary  pressure,  Pc,  is  defined  as  the  difference 
between  the  pressure  in  the  NAPL  phase  and  the  pressure  in  the  water  phase: 

P  =P  -P  [5] 

Darcy’s  law  can  be  written  for  each  phase  of  a  two-phase  immiscible  system,: 

,.=Z*k(vp,-p,gV2)  [6] 

kJc 

-  P„„gVz)  [7] 

where  w  denotes  water,  nw  denotes  NAPL,  q  is  the  volumetric  flux  (assumed  positive 
down),  p,  is  the  dynamic  viscosity,  k  is  the  intrinsic  permeability,  and  kr  is  the  relative 
permeability  for  each  phase. 

The  capillary  pressure  gradient  is  solved  for  by  combining  Equations  [5],  [6],  and 

[7]: 
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[8] 


VP  =  +  Ayojesina 


where  Ap=pnw  -  pw ,  and  a  is  the  angle  of  the  flow  from  the  horizontal.  The  product  of  k 
and  krw  is  kw  and  likewise  for  the  NAPL  phase. 

Equation  [8]  is  similar  to  Equation  [4].  The  factor,  sina,  accoimts  for  a  reduced 
gravity  effect  when  the  fluid-fluid  displacement  is  not  vertical.  This  term  can  likewise 
appear  in  Equation  [4].  The  important  difference  is  that  Equation  [4]  is  a  pore  scale 
equation  whereas  Equation  [8]  is  a  macroscopic  equation.  Equation  [4]  describes  the 
change  in  capillary  pressure  on  a  pore  scale.  The  change  in  capillary  pressure  on  the  pore 
scale  directly  indicates  the  stability  of  the  flow  in  each  pore  and  therefore  the  movement 
of  the  front.  Equation  [8]  has  the  same  form,  but  the  length  and  the  flux  are  defined  on  a 
macroscopic  scale,  and  k  is  a  macroscopic  permeability.  Equation  [8]  can  be  used  in  the 
determination  of  front  stability  (Hill,  1952;  Kueper  and  Frind,  1988),  but  using  a 
macroscopic  flux  in  an  equation  describing  microscopic  events  introduces  a  scaling  factor 
(Baudet  et  al.,  1986;  Oreskes  et  al.,  1994).  The  experimental  work  performed  by  Morrow 
and  Songkran  (1981)  and  Dawson  and  Roberts  (1997)  confirm  the  need  for  the  scaling 
factor,  however  neither  study  has  offered  an  explanation  for  its  derivation  or  magnitude. 
One  dimensional  air-NAPL  displacement  experiments  were  performed  by  Morrow  and 
Songkran  (1981)  to  determine  the  residual  saturation  of  the  nonwetting  fluid.  The 
residual  saturation  was  predicted  by  a  linear  combination  of  the  microscopic  capillary  and 
bond  numbers  when  the  capillary  number  was  increased  by  a  factor  of  1  O'*.  By 
increasing  the  capillary  number  by  a  factor  of  10'*,  the  capillary  and  bond  numbers  were 
of  the  same  order  of  magnitude  or  the  capillary  number  was  one  order  of  magnitude 
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higher  than  the  bond  number.  In  the  study  performed  by  Dawson  and  Roberts  (1997), 
one  dimensional  DNAPL-water  experiments  were  performed  and  the  DNAPL  saturation 
was  predicted  by  a  linear  combination  of  bond  and  capillary  numbers  when  the  capillary 
number  was  increased  by  the  inverse  of  the  relative  permeability  of  water.  In  the  case  of 
Dawson  and  Roberts  (1997),  the  use  of  the  relative  permeability  as  a  scaling  factor 
resulted  in  capillary  nmnbers  that  were  of  the  same  order  of  magnitude  as  the  bond 
number  or  one  order  of  magnitude  higher.  Both  studies  applied  the  linear  combination  of 
capillary  and  bond  numbers  to  an  examination  of  the  NAPL  residual  saturation.  Other 
researchers  (Wilson  et  al.,  1990)  have  used  the  linear  combination  in  addition  to  static 
force  balance  considerations  to  estimate  NAPL  mobility. 

The  stochastic  aggregation  algorithm  is  a  pore  scale  model  and  the  movement  of 
the  interface  is  described  by  Equation  [4].  The  problem  with  using  Equation  [4]  lies  in 
determining  the  values  of  pore  scale  parameters.  Instead  of  estimating  pore  scale 
parameters,  the  scale  incongruity  is  addressed  in  Equation  [8].  A  scaling  factor  of  10  was 
applied  to  the  wetting  and  nonwetting  capillary  numbers  for  the  current  study.  This  value 
was  chosen  because  it  yields  a  capillary  number  which  is  of  the  same  order  of  magnitude, 
or  one  order  of  magnitude  higher  than  the  bond  number  for  the  chosen 
nondimensionalization  following  die  results  of  Morrow  and  Songkran  (1981)  and 
Dawson  and  Roberts  (1997). 

The  length  scale  and  capillary  pressure  in  Equation  [8]  are  non-dimensionalized 
with  macroscopic  parameters  following  Hilfer  and  0ren  (1996).  A  hat  designates  the 
non-dimensional  variables: 

V  =  V/  [9] 
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[10] 


Non-dimensionalizing  the  capillary  pressure  gradient  results  in  Equation  [1 1]: 


^  QwMja  ,  A/^/sing  I 

kj*d  KJ*d  Pd 

where  a  accounts  for  the  change  from  pore  scale  to  the  macroscale  and  is  assigned  the 
value  of  10^  following  the  discussion  on  the  previous  page. 

The  first  two  terms  in  Equation  [1 1]  are  the  macroscopic  capillary  number  for 
water  and  NAPL,  respectively.  The  third  term  is  the  macroscopic  bond  number  as 
defined  by  Hilfer  and  0ren  (1996).  Equation  [1 1]  is  rewritten  in  terms  of  these 
dimensionless  numbers. 


VP,  =  Cofw  -  Canw  +  Bo  [12] 


where  Caw  and  Canw  are  the  capillary  numbers  for  water  and  NAPL,  respectively,  and  Bo 
is  the  bond  number.  The  bars  denote  that  the  dimensionless  numbers  are  macroscopic. 
Equation  [12]  quantifies  the  competition  between  the  capillary,  viscous,  and  gravity 
forces  at  the  interface. 

The  scaled  capillary  pressure  gradient  in  Equation  [12]  is  termed  the  transition 
number,  T,  and  is  grouped  as  follows; 

r  =  [13] 

The  first  term  of  the  transition  number  compares  the  macroscopic  capillary  numbers  for 
the  wetting  and  non-wetting  phases  and  captures  the  combined  effect  of  the  viscous  and 
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capillary  forces.  The  first  term  of  the  equation  is  positive  for  NAPLs  that  are  less  viscous 
than  water  and  negative  for  NAPLs  that  are  more  viscous  than  water.  The  second  term  in 
the  transition  number  is  the  macroscopic  bond  number  and  includes  the  effects  of  gravity 
forces  due  to  the  density  difference  between  the  NAPL  and  water.  The  second  term  is 
positive  for  DNAPLs  and  negative  for  LNAPLs,  if  the  NAPL  is  displacing  water  from 
the  top. 

3.1.2  Stochastic  Aggregation  Algorithm 

The  stochastic  aggregation  model  follows  a  simple  algorithm.  When  the 
simulation  starts,  particles  are  released  fi'om  a  random  location  opposite  the  seed  point 
(the  location  where  the  DNAPL  is  injected).  The  particle  has  an  equal  probability  of 
moving  up,  down,  right,  or  left.  If  the  particle  encounters  the  sides  of  the  network,  it 
“bounces  off’  and  continues  walking.  The  particle  continues  to  walk  randomly  until  it 
reaches  a  location  adjacent  to  a  location  already  occupied  by  a  particle.  The  particle  is 
incorporated  into  the  aggregate  with  a  probability  equal  to  the  sticking  probability,  and  a 
new  particle  is  released.  Once  incorporated  into  the  aggregate,  each  particle  represents 
the  DNAPL  occupying  a  pore  after  having  displaced  the  water  in  the  pore.  The  sticking 
probability  changes  the  thickness  of  the  branches  of  the  aggregate.  The  overall  density  of 
the  aggregate  is  increased  or  decreased.  This  effect  is  seen  by  comparing  the  two 
simulations  in  Figure  3.2.  The  sticking  probability  for  Figure  3.2(a)  is  67  percent  and  the 
sticking  probability  for  Figure  3.2(b)  is  2  percent.  The  change  in  aggregate  density  is  a 
function  of  the  grid  size  used. 
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3.2(a)  3.2(b) 

Figure  3.2.  Aggregates  formed  by  stochastic  aggregation  model  simulations  with  different 
sticking  probabilities.  The  sticking  probability  for  3.2(a)  and  3.2(b)  is  67%  and  2%,  respectively. 


The  model  is  calibrated  such  that  the  changes  in  the  front  stability  (viscous  fingers 
to  flat  fronts)  correspond  to  the  changes  in  the  aggregate  density  as  the  sticking 
probability  is  increased  or  decreased.  The  transition  number  determines  the  changes  in 
front  stability.  The  transition  number  is  related  to  the  sticking  probability  via  calibration 
with  laboratory  experiments. 


3.2  Model  Calibration 

A  series  of  two-dimensional  laboratory  experiments  involving  DNAPLs 
displacing  water  in  a  glass  bead  porous  medium  were  performed  at  Colorado  State 
University.  The  flumes  were  approximately  36  cm  in  width,  35  cm  in  length,  and  1  glass 
bead  in  thickness.  A  DNAPL  was  injected  into  the  water-saturated  media  through  a  tube 
at  the  top  of  the  flume  at  a  constant  flow  rate.  The  flow  rate  and  the  angle  of  inclination 
of  the  flvime  from  the  horizontal  were  systematically  varied.  The  flux  of  the  DNAPL 
entering  the  soil  was  determined  based  on  the  reported  flow  rate  and  the  area  of  flow  of 
the  tubing  used  for  injection.  The  cross  sectional  area  of  the  flume  was  not  used  to 
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determine  the  flux  for  the  calibration,  and  neither  the  flume  width  nor  the  tubing  width 
provides  the  correct  dimension  for  the  flux  determination.  The  flume  width  is  an 
arbitrary  dimension  that  is  unknown  in  a  field  application  of  the  model;  however,  the  area 
through  which  the  DNAPL  was  leaked  may  be  estimated.  Each  experiment  was 
photographed  at  different  times  starting  fi-om  the  time  of  the  DNAPL  injection.  Carbon 
Tetrachloride  (CTET),  1,2  Dichloroethane  (DC  A),  Tetrachloroethylene  (PERC),  and 
Mobile  P5^ogard  53  (made  by  Akzo  Nobel,  Gallipolis  Ferry,  WV,  vmder  the  name 
Fyrquel  220)  were  the  DNAPLs  used  in  the  experiments.  The  first  three  are  commonly 
used  DNAPLs  which  are  less  viscous  than  water.  Mobile  Pyrogard  53  is  a  fire-resistant 
hydraulic  fluid  and  is  more  viscous  than  water. 

Collectively,  the  laboratory  experiments  show  that: 

•  for  a  DNAPL  which  is  less  viscous  than  water,  the  front  becomes  more 
unstable  as  both  the  NAPL  flow  rate  and  the  angle  from  the  horizontal  are 
increased. 

•  for  a  DNAPL  which  is  more  viscous  than  water,  the  firont  becomes  more  stable 
as  the  NAPL  flow  rate  increases,  but  the  fi'ont  becomes  less  stable  as  the  angle 
from  the  horizontal  increases. 

The  experimental  photographs  were  digitized  and  analyzed  to  determine  the  ratio 
of  DNAPL-occupied  area  to  total  area  using  GLOBAL  LAB®  Image  (Copyright  ©  1995 
by  Data  Translation,  Inc.,  Marlboro,  MA),  a  Windows-based  image  processing  and 
analysis  package.  Forty-three  laboratory  experiments  were  used  to  calibrate  the  model. 
The  calibration  group  included  three  different  volumetric  flow  rates  (9ml/min,  15ml/min, 
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and  30  ml/min),  where  each  DNAPL  was  represented  at  different  angles  of  inclination. 

The  parameters  and  pictures  for  each  model  calibration  case  are  available  in  Appendix  B. 

The  stochastic  aggregation  model  was  run  to  determine  the  average  ratio  of 
occupied  area  to  total  area  for  each  sticking  probability.  The  model  was  stopped  when 
the  particles  reached  the  bottom  of  the  grid.  The  grid  size  was  calculated  to  be  1 13  X 
135,  where  each  grid  space  represented  1  pore  in  the  experimental  photographs.  The 
calibration  of  the  stochastic  aggregation  model  is  unique  to  this  grid  size.  The  grid  used 
for  the  simulations  must  be  square  in  order  to  use  the  graphics  program,  PSPLOT 
(written  by  Kevin  Kohler,  Nova  Southeastern  University  Oceanographic  Center,  Dania, 
FL)  to  transform  the  model  output.  This  is  a  limitation  in  PSPLOT  and  not  with  the 
model  itself  Because  of  the  limitation  in  the  graphics  program,  all3X113  grid  was 
used.  The  vertical  proportions  were  kept  the  same  as  the  laboratory  simulations  but  the 
sides  were  truncated.  The  only  laboratory  experiments  where  the  DNAPL  approaches  the 
sides  of  the  flume  are  those  involving  Mobile  Pyrogard  53.  The  model  was  run  16  times 
each  for  23  sticking  probabilities  (ranging  from  Ps=0.0001  to  Ps=l)  to  determine  an  area 
ratio.  The  area  ratio  is  defined  as  the  average  ratio  of  occupied  area  to  total  area  for  each 
sticking  probability.  One  simulation  runs  in  0.5  to  15  minutes  (real  time)  on  an  IBM 
RS6000  depending  on  the  value  of  the  sticking  probability.  The  area  ratio  is  a  function  of 
the  transition  number  for  the  laboratory  experiments.  For  the  stochastic  aggregation 
model  the  area  ratio  is  a  function  of  the  sticking  probability  for  the  selected  grid  size. 

The  relationship  between  the  sticking  probability  and  the  area  ratio  for  the 
stochastic  aggregation  model  is  shown  in  Figure  3.3.  The  data  is  fit  with  the  following 
power  law  relationship: 
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model  calibration  with  the  laboratory  experiments  was  0.376,  and  is  on  the  portion  of  the 
curve  where  a  high  level  of  the  area  ratio’s  variability  is  accounted  for  by  the  sticking 
probability. 

Equation  [14]  is  rearranged  so  that  the  area  ratio  (AR)  is  the  independent  variable. 

Ps  =  1 .04323  X  10-^(y4i?)-^ 

The  AR  determined  from  the  photographs  of  laboratory  experiments  was  used  as  input  to 
Equation  [15]  to  determine  a  sticking  probability  for  each  photograph.  A  relationship 
between  the  transition  number  and  the  sticking  probability  is  then  determined  by  plotting 
the  values  for  each  of  the  lab  experiments  in  the  chosen  calibration  group.  Figure  3.4 
shows  the  relationship  between  the  absolute  value  of  the  transition  number  and  the 
sticking  probability  as  well  as  the  corresponding  power  law  fit.  Seven  of  the  laboratory 
experiments  had  area  ratios  that  were  low  enough  to  yield  a  sticking  probability  of  greater 
than  one  using  Equation  [15].  These  data  points  are  shown  on  Figure  3.4,  but  the  power 
law  relationship  is  dotted  in  this  region  as  the  model  cannot  use  a  sticking  probability 
greater  than  one  as  input.  Using  the  absolute  value  of  the  transition  number  allows  the 
model  to  use  small  sticking  probabilities  to  capture  phenomena  resulting  from  two 
different  mechanisms.  One  mechanism  is  the  change  in  the  DNAPL-water  displacement 
resulting  from  changes  in  intrinsic  permeability  within  the  porous  media.  The  intrinsic 
permeability  (k)  may  be  varied  in  the  calculation  of  the  transition  number  to  account  for 
porous  media  heterogenities.  Variations  in  k  are  included  in  the  comparisons  with 
heterogeneous  laboratory  experiments  in  Chapter  5.  In  this  case,  a  low  sticking 
probability  prohibits  the  particles  from  entering  low  permeability  regions  of  the  grid. 
Large  values  for  the  transition  number  are  seen  when  the  intrinsic  permeability  of  the 
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The  sticking  probability  was  calculated  for  the  laboratory  experiments  not 
included  in  the  calibration  group  using  Equation  [16].  Model  simulations  corresponding 
to  each  of  these  experiments  are  discussed  in  greater  detail  in  Chapter  4. 
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CHAPTER  4 


MODEL  CONFIRMATION  AND  SCALING  IN  HOMOGENEOUS 

POROUS  MEDIA 


Models  that  attempt  to  predict  natural  phenomena  are  used  in  every  area  of 
science  and  engineering.  The  challenge  inherent  in  developing  a  model  to  represent  real- 
life  processes  is  in  assessing  the  model’s  adequacy.  The  terms  verification  and  validation 
are  often  used  to  make  claims  of  a  model’s  adequacy.  Although  techniques  for  verifying 
and  validating  models  have  been  developed,  it  is  impossible  to  truly  verify  or  validate  a 
numerical  model  for  a  natural  system  (Oreskes  et  al.,  1994).  Verification  of  a  model  is 
only  possible  in  a  closed  system  in  which  all  system  components  are  established 
independently.  Validation  of  a  model  means  that  the  model  is  legitimate.  All  input 
parameters  as  well  as  the  assumptions  used  must  be  known  to  be  accurate  to  show  that  a 
model  has  been  validated. 

Even  though  models  cannot  be  verified  or  validated,  their  adequacy  must  be 
assessed  in  some  way.  The  best  we  can  do  to  assess  a  model’s  adequacy  is  to  compare 
the  model  results  with  laboratory  or  field  data  and  call  this  process  model  confirmation. 
In  model  confirmation,  laboratory  or  field  observations  which  match  model  output  are 
said  to  support  the  probability  that  the  model  has  been  formulated  correctly  (Oreskes  et 
al.,  1994). 
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4.1  Model  Confirmation 

Progress  towards  model  confirmation  is  accomplished  by  showing  that  model 
realizations  match  laboratory  or  field  experiments  for  DNAPL-water  displacement.  In 
this  chapter,  the  stochastic  aggregation  model  simulations  are  compared  to  two- 
dimensional  laboratory  DNAPL-water  experiments  in  homogeneous  porous  media.  In 
Chapter  5,  model  simulations  are  compared  with  heterogeneous  systems. 

Model  comparisons  consist  of  33  laboratory  DNAPL  experiments.  The  DNAPLs 
in  the  validation  group  are  the  same  as  those  in  the  calibration  group,  however  the 
volumetric  flow  rates  are  different.  Table  4.1  lists  the  DNAPL,  flow  rate,  flux,  time  at 
which  the  picture  was  taken,  transition  number,  corresponding  sticking  probability,  area 
ratio  as  calculated  from  the  photographs  (ARi),  and  the  area  ratio  from  the  model 
simulations  (AR2).  The  area  ratio  for  the  model  simulations  is  an  average  area  ratio  to 
within  one  half  of  the  standard  deviation  using  a  95%  confidence  interval  and  16  model 
runs.  In  six  of  these  cases  the  transition  number  corresponds  to  a  sticking  probability  of 
greater  than  one.  The  sticking  probability  is  greater  than  one  when  the  transition  number 
is  less  than  or  equal  to  0.205 141 .  This  occurs  when  a  combination  of  a  low  flow  rate  and 
a  small  angle  of  inclination  from  the  horizontal  yields  a  low  value  for  the  transition 
number.  A  sticking  probability  of  one  is  used  in  these  cases,  as  an  approximation,  with 
good  qualitative  results.  The  NAPL  did  not  reach  the  bottom  of  the  flume  in  the  Mobile 
Pyrogard  53  laboratory  experiments.  The  volume  of  DNAPL  released  is  determined  to 
make  possible  a  comparison  between  the  lab  experiments  and  model  simulations.  The 
volume  is  related  to  a  model  simulation  volume  or  number  of  particles  based  on  the 
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approximate  pore  size.  The  model  is  stopped  when  the  calculated  number  of  particles  is 
incorporated  into  the  simulated  aggregate. 


Table  4.1 

Model  Confirmation  Cases 


DNAPL 

Angle 

Q 

q 

Time 

T 

Ps 

ARi  ARj 

■fflirniiiifl 

(min) 

1,2  DCA 

90“ 

3 

0.0617 

1 

0.393 

0.423 

0.0424  0.0873 

1,2  DCA 

90“ 

6 

0.1234 

1 

0.544 

0.275 

0.0585  0.1040 

1,2  DCA 

45“ 

3 

0.0617 

1 

0.322 

0.550 

0.0731  0.0785 

1,2  DCA 

45“ 

6 

0.1234 

1 

0.473 

0.331 

0.145  0.0965 

1,2  DCA 

15“ 

3 

0.0617 

2 

0.214 

0.948 

0.0804  0.0629 

1,2  DCA 

15“ 

6 

0.1234 

1 

0.364 

0.467 

0.1085  0.0838 

CTET 

90“ 

3 

0.0617 

1 

0.313 

0.571 

0.0644  0.0773 

CTET 

90° 

6 

0.1234 

1 

0.338 

0.516 

0.1003  0.0805 

CTET 

45° 

3 

0.0617 

1 

0.229 

0.866 

0.0612  0.0653 

CTET 

45° 

6 

0.1234 

1 

0.254 

0.755 

0.0986  0.0669 

CTET 

30° 

3 

0.0617 

2 

0.169 

1.292' 

0.0445  0.0640 

CTET 

30° 

6 

0.1234 

1 

0.194 

1.078' 

0.0623  0.0640 

CTET 

15“ 

3 

0.0617 

1 

0.099 

2.608' 

0.0275  0.0640 

CTET 

15“ 

6 

0.1234 

1 

0.124 

1.942' 

0.0464  0.0640 

PERC 

90“ 

3 

0.0617 

1 

0.354 

0.485 

0.0551  0.0826 

PERC 

90“ 

6 

0.1234 

1 

0.402 

0.410 

0.0825  0.0883 

PERC 

90“ 

12 

0.2467 

0.5 

0.496 

0.310 

0.0913  0.0990 

PERC 

45“ 

3 

0.0617 

1 

0.264 

0.715 

0.0645  0.0706 

PERC 

45“ 

6 

0.1234 

1 

0.312 

0.574 

0.0636  0.0771 

PERC 

45“ 

12 

0.2467 

1 

0.407 

0.404 

0.1172  0.0889 

PERC 

30“ 

3 

0.0617 

1 

0.207 

0.991 

0.0402  0.0648 

PERC 

o 

o 

6 

0.1234 

1 

0.248 

0.777 

0.0344  0.0682 

PERC 

15“ 

3 

0.0617 

1 

0.127 

1.890' 

0.0295  0.0640 

PERC 

15“ 

6 

0.1234 

1 

0.174 

1.241' 

0.0429  0.0640 

MP53^ 

90“ 

3 

0.0309 

0.08 

-2.44 

0.0377 

0.0060  0.0056 

MP53^ 

90“ 

3 

0.0309 

0.5 

-2.44 

0.0377 

0.0215  0.0332 

MP53^ 

90° 

3 

0.0309 

0.75 

-2.44 

0.0377 

0.0311  0.0498 

MP53^ 

90“ 

3 

0.0309 

1 

-2.44 

0.0377 

0.0367  0.0665 

MP53^ 

90“ 

6 

0.0619 

1 

-4.98 

0.0146 

0.0742  0.1329 

MP53^ 

90“ 

6 

0.0619 

3 

-4.98 

0.0146 

0.2137  0.3987 

MP53^ 

45“ 

3 

0.0309 

1 

-2.47 

0.0371 

0.0708  0.0665 

MP53^ 

45“ 

6 

0.0619 

1 

-5.02 

0.0145 

0.0664  0.1329 

MP53^ 

45“ 

6 

0.0619 

3 

-5.02 

0.0145 

0.1734  0.3987 

1 .  A  sticking  probability  of  1  was  used  in  these  cases. 

2.  In  these  experiments  the  DNAPL  did  not  reach  the  bottom  of  the  flume,  and  the  model  was  adjusted 
for  the  corresponding  volume. 
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The  area  ratio  (ARi)  calculated  from  the  photographs  of  the  laboratory 
experiments  is  plotted  on  the  x-axis  and  the  average  area  ratio  (AR2)  from  18  model 
simulations  is  plotted  on  the  y-axis  in  Figure  4.1.  Figure  4.1  shows  the  model’s 
performance  in  producing  simulations  wherein  the  amount  of  DNAPL  in  the  system  at  a 
given  flow  rate  and  time  may  be  compared  to  the  laboratory  experiments  using  the  value 
of  the  area  ratio.  The  data  points  for  Mobile  Pyrogard  53  are  excluded  from  Figure  4.1 
because  their  model  simulated  area  ratios  (AR2)  were  pre-determined.  The  error  bars 
represent  a  90%  confidence  interval  for  the  model  simulations.  The  data  points  would  lie 
on  the  1 : 1  line  shown  on  Figure  4. 1  if  the  model  simulations  matched  the  experimental 
photographs  perfectly.  There  is  no  clear  trend  indicating  that  the  model  consistently 
over-  or  under-predicts  the  area  ratio  of  the  laboratory  experiments.  The  variability  can 
be  attributed  to  the  error  in  the  conditions  under  which  the  laboratory  experiments  were 
performed,  as  well  as  error  in  evaluating  the  area  ratio  based  on  the  photographs.  Using 
Figure  4.1  as  the  only  means  of  evaluating  the  model  assumes  that  both  the  experiments 
and  the  area  ratio  determination  from  the  photographs  are  accurate.  Figure  4.1  does  not 
give  any  information  regarding  the  characteristic  way  in  which  the  DNAPL  displaces  the 
water.  The  characteristic  way  in  which  a  DNAPL  displaces  water  has  been  a  primary 
focus  of  this  research  and  is  evaluated  by  visual  comparison. 
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Area  Ratio  Lab  Experiments  (AR1) 

Figure  4.1.  A  comparison  of  the  area  ratios  as  calculated  from  the  lab  experiments 

and  the  model  simulations. 

Two  data  points  furthest  from  the  1:1  line  belong  to  DCA  simulations.  For  DCA 
at  90°,  6  ml  min  \  and  1  minute  (point  [0.0585,  .01040])  the  model  predicts  a  higher  area 
ratio  than  was  observed  in  the  laboratory  experiment.  For  DCA  at  45°,  6  ml  min"*,  and  1 
minute  (point  [0.145, 0.0965])  the  model  predicts  a  lower  area  ratio  than  was  observed  in 
the  laboratory.  The  inconsistency  is  attributed  to  the  poor  contrast  of  the  photographs 
used  in  determining  the  area  ratio  (see  Appendix  C).  When  the  contrast  of  the 
photograph  is  poor,  the  software  cannot  easily  differentiate  between  areas  where  the 
DNAPL  is  present  and  darkened  areas  of  the  surrounding  porous  media. 

The  effects  of  increasing  the  flow  rate  and  increasing  the  angle  of  inclination  are 
examined  for  the  DNAPLs  used  in  the  model  confirmation.  Increasing  the  angle  of 
inclination  from  the  horizontal  and  increasing  the  flow  rate  both  increase  the  instability  of 
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the  DNAPL-water  displacement  if  the  DNAPL  is  less  viscous  than  water.  In  the  case  of  a 
DNAPL  which  is  more  viscous  than  water,  increasing  the  flow  rate  increases  the  stability 
of  the  DNAPL-water  front  while  an  increase  in  the  angle  of  inclination  increases  the 
instability. 

Smaller  values  of  the  transition  number  indicate  more  stable  fronts  for  the  range 
of  transition  numbers  in  Table  4.1 .  This  is  seen  in  the  case  of  PERC  for  two  different 
angles  of  inclination  and  a  flow  rate  of  3  ml  min'^  At  90°,  the  value  of  the  transition 
number  is  0.354.  With  a  decrease  in  the  angle  of  inclination  from  the  horizontal  to  15°, 
the  transition  number  reduces  to  0.127.  The  same  trend  is  seen  if  the  flow  rate  is 
changed.  For  PERC  with  an  angle  of  inclination  of  90°,  the  transition  number  increases 
from  0.354  to  0.402  as  flow  rate  increases  from  3  ml  min'^  to  6  ml  min'. 

Whereas  DCA,  CTET,  and  PERC  are  less  viscous  than  water.  Mobile  Pyrogard 
53  is  more  viscous  than  water.  Because  Mobile  Pyrogard  53  is  more  viscous  than  water, 
the  difference  in  the  viscosities  used  to  calculate  the  transition  number  becomes  negative. 
The  influence  of  the  flux  rate  on  the  magnitude  and  sign  of  the  transition  number  is 
shown  by  changing  the  flow  rate  from  3  ml  min'*  to  6  ml  min'*  for  Mobile  Pyrogard  53  at 
90°.  The  transition  number  decreases  from  —2.44  to  —4.98  with  the  increase  in  flux.  The 
magnitude  of  T  for  Mobile  Pyrogard  53  is  large  because  of  the  absolute  value  of  the 
viscosity  difference  of  4.86  cP  between  Mobile  Pyrogard  53  and  water.  The  viscosity 
difference  with  water  ranges  from  0.1 13  cP  to  0.035  cP  for  the  other  DNAPLs  used  in  the 
experiments.  The  increase  in  the  angle  of  inclination  from  the  horizontal  results  in  larger 
values  of  the  transition  number  for  Mobile  Pyrogard  53.  This  is  seen  by  looking  at 
Mobile  Pyrogard  53  at  a  flux  of  3  ml  min'*.  The  transition  number  increases  from  -2.47  to 


30 


-2.44  for  an  increase  in  the  angle  of  45®  to  90®  from  the  horizontal.  This  increase  is 
small,  but  follows  the  general  trend.  Tables  1  through  4  in  Appendix  C  contain  the  input 
values  required  for  the  calculation  of  the  transition  number  for  each  case  listed  in  Table 
4.1. 

The  sticking  probability  approaches  positive  infinity^  positive  transition 
numbers  approach  zero.  This  leads  to  a  definition  of  the  model’s  limits.  The  sticking 
probability  is  greater  than  one  when  the  transition  number  is  less  than  or  equal  to 
0.205 141 .  Combinations  of  inputs  resulting  in  the  transition  number  being  less  than 
0.205 141  form  a  limit  for  the  model’s  applicability.  This  limit  is  encountered  in  the 
model  validation  cases  PERC  and  CTET  with  of  a  low  angle  of  inclination  and  a  low 
flux.  The  area  ratio  determined  by  the  laboratory  experiment  is  smaller  than  the  smallest 
area  ratio  that  the  model  could  provide  at  a  sticking  probability  of  1 .  As  the  absolute 
value  of  the  transition  number  goes  to  infinity,  the  sticking  probability  approaches  zero. 
The  limit  of  the  model  as  the  transition  number  increases  depends  on  the  amount  of  time 
available  for  a  simulation.  This  is  determined  by  the  computer  or  by  physical 
circumstances. 

Figures  4.2  through  4.6  show  examples  of  comparisons  between  the  laboratory 
experiments  and  the  model  simulations.  A  complete  set  of  model  validation  pictures  is 
shown  in  Appendix  C.  Four  sets  of  pictures  are  shown  for  DCA,  CTET  and  PERC 
including  pictures  at  two  different  angles  of  inclination  and  two  different  inflow  rates. 
Eight  pictures  are  shown  for  Mobile  Pyrogard  53.  Four  of  the  Mobile  Pyrogard  53 
pictures  show  the  model  simulations  as  time  progresses  for  the  same  angle  of  inclination 
and  inflow  rate.  The  other  four  Mobile  Pyrogard  53  pictures  show  changes  in  the  front 
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with  changes  in  the  angle  of  inclination  and  inflow  rate.  It  is  important  to  keep  in  mind 
that  the  objective  of  this  type  of  modeling  is  not  an  exact  picture  of  how  the  DNAPL 
displaced  the  water,  but  a  characteristic  picture  of  how  the  DNAPL  displaced  the  water. 
Characteristics  such  as  the  thickness  of  a  finger,  the  number  of  fingers  and  the  relative 
amount  of  pore  space  occupied  are  important  to  compare. 

The  model  simulations  and  experimental  pictures  match  quite  well.  In  the  cases 
of  DCA,  CTET  and  PERC  (Figures  4.2, 4.3,  and  4.4),  the  branching  and  thickness  of  the 
fingers  are  comparable.  At  90°  and  3  ml  min  ‘  the  respective  DNAPL  appears  to  go 
straight  down  through  the  porous  media  with  little  branching.  When  the  flow  rate  is 
doubled  to  6  ml  min  \  more  branching  can  be  seen  as  the  DNAPL  travels  through  the 
flume.  As  the  flow  rate  increases,  the  front  is  increasingly  unstable  for  DCA,  CTET,  and 
PERC.  Often,  increased  instability  is  thought  to  coincide  with  thiimer,  deeper  fingers. 
Both  the  experimental  photographs  and  the  model  simulations  show  that  this  is  not  true. 
The  capillary  pressme  gradient  increases  as  the  DNAPL  travels  deeper,  leading  to  the 
fi'ont  becoming  increasingly  imstable.  If  the  capillary  pressure  gradient  increases,  then 
the  capillary  pressure  at  the  fi'ont  must  increase  eis  the  distance  to  the  front  increases. 
More  pores  are  invaded  leading  to  more  lateral  movement  at  higher  flow  rates  as  the 
capillary  pressure  at  the  firont  increases. 

In  Figure  4.3  for  CTET  at  15°,  the  model  simulations  in  (c)  and  (d)  have  a 
sticking  probability  of  one.  In  both  of  these  cases,  the  calculated  sticking  probabilities 
were  greater  than  one.  A  sticking  probability  of  one  produced  model  realizations  which 
compare  well  with  the  laboratory  experiment. 
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The  comparison  between  the  simulations  and  photographs  for  Mobile  Pyrogard  53 
are  likewise  very  good.  Figure  4.5  shows  Mobile  Pyrogard  53  displacing  water  at  3  ml 
min*  at  5  seconds,  30  seconds,  45  seconds,  and  1  minute.  The  sticking  probability 
captures  the  characteristic  displacement  of  the  front.  The  model  is  stopped  when  the 
approximated  number  of  pores  are  filled  based  on  the  volume  determination.  This  does 
not  imply  diat  die  model  can  account  for  time  variation,  as  there  is  no  time  parameter  in 
the  equation  for  the  transition  number.  Other  generalized  DLA  models  have  incorporated 
a  temporal  element  into  the  algorithm  (Flury  and  Fluhler,  1995).  The  current  model  can 
show  realizations  at  different  times  if  the  volume  is  known.  Figure  4.6  shows  the  effect 
of  changing  the  angle  of  inclination  from  90°  to  45°  and  the  flow  rate  from  3  ml  min  *  to 
6  ml  min  *  at  1  minute.  The  Mobile  Pyrogard  53  seems  to  travel  farther  in  the  simulation 
than  in  the  experiment  in  Figures  4.6(b)  and  4.6(d).  This  is  due  to  the  approximation  of 
the  number  of  pores  and  the  pore  size.  The  same  size  grid  is  used  for  all  of  the  model 
simulations,  but  the  experiments  were  performed  by  repacking  the  flume  for  each 
scenario.  It  is  difficult  to  determine  the  nature  of  the  packing  for  each  experiment  after 
the  fact.  The  number  of  particles  corresponding  to  a  particular  volume  is  approximated 
based  on  an  average  pore  size.  Both  of  these  facts  contribute  to  uncertainty  in  calculating 
the  exact  number  of  particles  corresponding  to  a  given  volume.  This  issue  will  be 
brought  up  again  in  Chapter  5  in  regard  to  natural  sands.  The  transition  numbers  for 
Mobile  Pyrogard  53  indicate  that  the  DNAPL-water  front  is  stable,  as  the  capillary 
pressure  gradient  is  negative.  The  sticking  probabilities  corresponding  to  the  negative 
transition  number  are  small  leading  to  model  simulations  with  more  compact  fronts. 
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(a)  90°,  0=3  ml  min  ,  time=l  min,  T=0.393,  Ps=0.422 


(b)  90°,  0=6  ml  min',  time=l  min,  T=0.544,  Ps=0.275 
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(c)  15°,  0=3  ml  min'',  time=2  min,  T=0.214,  Ps=0.948 


(d)  15°,  0=6  ml  min*',  time=l  min,  T=0.364,  Ps=0.467 


Figure  4.2.  Lab  experiments  and  model  simulations  for  1,2  Dichloroethane  (DCA). 


34 


(a)  90®,  Q=3  ml  min"',  time=l  min,  T=0.313,  Ps=0.571 


(c)  15®,  0=3  ml  min'*,  time=l  min,  T=0.099,  Ps=2.608' 
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(d)  15®,  0=6  ml  min'*,  time=l  min,  T=0.124,  Ps=1.942' 


Figure  4.3.  Lab  experiments  and  model  simulations  for  Carbon  Tetrachloride  (CTET). 
‘Note:  A  sticking  probability  of  1.0  was  used  in  the  model  simulation. 
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(a)  90®,  0=3  ml  min'*,  time=l  min,  T=0.354,  Ps=0.485 


(a)  90®,  Q=3  ml  min'*,  time=5  s,  T=-2.44j  Ps=0.0377 


(b)  90°,  Q=3  ml  min'*,  time=30  s,  T—2.44,  Ps=0.0377 


Figure  4.5.  Lab  experiments  and  model  simulations  for  Mobile  Pyrogard  53  at  different  times. 
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(a)  90°,  0=3  ml  min'*,  time=l  min,  T=-2.44,  Ps=0.0377 


(b)  90°,  Q=6  ml  min'',  time=l  min,  T=-4.98,  Ps=0.0146 


(c)  45°,  Q=3  ml  min"',  time=l  min,  T=-2.47,  Ps=0.0371 


(d)  45°,  Q=6  ml  min'',  time=l  min,  T=-5.02,  Ps=0.0145 


Figure  4.6.  Lab  experiments  and  model  simulations  for  Mobile  Pyrogard  53. 
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4.2  Model  Scaling 

Model  confirmation  was  performed  using  a  grid  size  corresponding  to  the 
approximate  number  of  pores  in  the  flume.  The  changes  in  the  displacement  fronts 
resulting  from  changes  in  flow  rate  and  angle  of  inclination  from  the  horizontal  are 
apparent  even  though  these  flumes  were  small.  Other  laboratory  experiments  found  in 
the  literature  (Illangasekare  et  al.,  1995  and  Kueper  and  Frind,  1991)  used  larger  flumes, 
and  field  spills  are  at  an  even  larger  scale.  It  is  of  interest  to  investigate  the  scalability  of 
the  model. 

Grids  up  to  an  order  of  magnitude  larger  than  the  model  validation  simulations  are 
investigated.  The  area  ratio  for  the  model  simulations  is  defined  as  the  number  of 
occupied  grid  spaces  divided  by  the  total  number  of  grid  spaces.  The  area  ratio  decreases 
with  increasing  grid  size  for  the  same  sticking  probability.  Figure  4.7  shows  the 
relationship  between  the  area  ratio  and  the  grid  size  for  a  sticking  probability  of  0.1.  The 
data  result  from  running  the  model  at  each  grid  size  for  the  given  sticking  probability. 

The  data  points  do  not  reflect  an  average  at  the  corresponding  grid  size.  The  time  it  takes 
to  complete  a  simulation  increases  dramatically  as  the  grid  size  increases.  The 
simulations  generated  for  model  validation  onall3xll3  grid  required  a  few  seconds  to 
15  minutes  (real  time)  depending  on  the  sticking  probability  and  the  user  availability  on 
an  IBM  RS6000.  The  simulation  time  increased  from  3  minutes  for  a  200  x  200  to  87 
minutes  for  a  500  x  500  with  a  sticking  probability  of  0.1 .  The  simulation  using  an  800  x 
800  grid  took  approximately  4  days.  It  is  prohibitive  to  run  multiple  simulations  and 
produce  average  area  ratios. 
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Figure  4.7.  The  decrease  in  area  ratio  with  increasing  grid  size  shown  for  a  sticking  probability  of  0.1. 

The  change  in  area  ratio  with  change  in  grid  size  does  not  prohibit  the  scaling  of 
the  model.  The  characteristics  of  the  aggregate  do  not  change  for  a  given  sticking 
probability  even  though  the  relative  amount  of  occupied  grid  space  is  changing.  This 
idea  is  similar  to  a  situation  in  the  field  where  a  DNAPL  displaces  water  in  a  narrow 
finger-like  pattern.  The  characteristic  of  the  finger  (thickness,  branching,  etc.)  would  be 
the  same  whether  the  volume  of  soil  immediately  surrounding  the  finger  is  considered  or 
a  much  larger  volume  is  considered.  The  ratio  of  occupied  pore  space  to  total  pore  space 
changes  as  the  selected  volume  changes.  Figure  4.8  shows  two  model  simulations  with 
the  same  sticking  probability  (Ps=0.1)  on  different  grid  sizes.  The  grid  size  for  4.8(a)  is 
200  X  200  and  the  grid  size  for  4.8(b)  is  600  x  600.  Simulations  for  grid  sizes  800  x  800 
and  1000  x  1000  were  run,  but  the  output  cannot  be  visualized  at  this  time.  The  darkness 
of  4.8(b)  is  due  to  a  limitation  in  the  program  used  to  generate  the  picture  from  the 
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model’s  output,  PSPLOT.  The  program  generates  the  picture  in  a  space  of  4.7  cm  by  4.7 
cm  regardless  of  the  model’s  grid  size.  The  pictures  become  darker  and  darker  with  the 
increasing  number  of  lines  being  drawn  in  the  same  space. 


4.8(a) 


4.8(b) 


Figure  4.8.  Model  simulations  for  a  sticking  probability  of  0.1.  Figure  4.8(a)  shows  a  model  simulation  on 
a  200  X  200  grid.  Figure  4.8(b)  shows  a  model  simulation  on  a  600  x  600  grid. 


The  question  of  model  scaling  is  addressed  by  determining  if  the  area  ratio  for  a 
given  sticking  probability  remains  characteristically  die  same.  This  implies  that  Figure 
4.8(a)  has  the  same  eirea  ratio  as  4.8(b)  if  a  200  x  200  portion  of  4.8(b)  is  examined.  A 
smaller  sized  grid  can  be  used  for  simulations  to  show  the  characteristic  DNAPL-water 
displacement  if  this  is  the  case.  Model  simulations  with  grid  sizes  ranging  over  one  order 
of  magnitude  are  examined.  The  results  are  presented  in  Table  4.2  and  Figure  4.9.  Np  is 
the  number  of  particles  occupying  spaces  inall3xll3  grid  centered  about  the  injection 
point  at  the  top  of  the  grid  (the  seed  particle).  The  simulations  are  stopped  when  the 
longest  finger  reached  row  113.  The  area  ratio  from  the  larger  grid  is  directly  compared 
to  the  average  area  ratio  forall3xll3  grid  at  the  same  sticking  probability  (Ps=0.01). 
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Table  4.2 

Area  ratio  for  a  1 13  x  1 13  sub-grid  on  increasing  grid  sizes  for  a 
sticking  probability  of  0.0 1 . 


Grid  Size 

Np 

Area  Ratio 

113x113 

1854 

0.1452 

200  X  200 

1939 

0.1518 

300  X  300 

1694 

0.1327 

400x400 

2339 

0.1832 

500  X  500 

1924 

0.1507 

800  X  800 

2026 

0.1587 

1000  X 1000 

2240 

0.1698 

Figure  4.9  shows  the  area  ratio  oscillating  and  perhaps  increasing  slightly  over  the 
range  of  grid  sizes  investigated.  One  possible  explanation  for  this  occxirrence  is  the 
increase  in  lateral  grid  space.  The  aggregate  ends  up  with  more  occupied  grid  spaces  in  a 
given  row,  increasing  the  area  ratio  slightly  when  there  is  more  lateral  space  in  which  the 
particle  can  move.  This  increase  is  small  and  results  in  a  1 .35%  increase  in  occupied 
pore  space  over  an  order  of  magnitude  change  in  grid  size.  It  is  concluded  from  this 
investigation  that  the  grid  may  be  scaled  to  the  113x113  size  used  in  the  model 
validation  with  little  change  to  the  characteristics  of  the  DNAPL-water  displacement  in 
the  pore  space.  This  conclusion  should  be  confirmed  by  studying  the  grid  size  and 
change  in  area  ratio  over  several  orders  of  magnitude  in  future  studies. 
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Grid  Size 


Figure  4.9.  The  area  ratio  for  a  1 13  x  1 13  portion  of  the  simulated  grid  compared 
for  grid  varying  sizes  up  to  1000  x  1000. 


The  scaling  technique  is  used  to  evaluate  the  basic  characteristics  of  larger 
laboratory  experiments  and  field  spills  by  scaling  the  DNAPL  volume.  The  volume  of 
Mobile  Pyrogard  53  is  related  to  the  number  of  particles  in  the  model  simulation  in  the 
model  validation  section.  This  is  accomplished  by  estimating  the  percentage  of  DNAPL 
occupied  pore  space  using  the  volume  of  DNAPL  injected  and  the  approximate  volume 
of  the  pore  space.  The  number  of  particles  necessary  for  a  simulation  of  a  given  volume 
is  determined  by  multiplying  the  percentage  of  occupied  pore  space  by  the  number  of 
pores  or  grid  spaces  for  the  simulation  (1 13^).  This  technique  is  utilized  in  the  simulation 
of  the  laboratory  experiments  of  Illangasekare  et  al.  (1995)  in  Chapter  5. 
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CHAPTERS 


MODEL  CONFIRMATION  FOR  HETEROGENEOUS  POROUS  MEDIA 

5.1  Heterogeneous  Model  Confirmation 

Heterogenities  are  incorporated  into  the  model  by  changing  the  transition  number. 
In  Chapter  4,  changes  in  the  transition  number  corresponding  to  changes  in  the  sticking 
probability  are  explored  Jis  a  function  of  the  fluid  properties  (density  and  viscosity)  of 
the  DNAPL  and  the  angle  of  inclination  from  the  horizontal.  The  characteristic  DNAPL- 
water  front  due  to  porous  media  heterogenieties  is  explored  by  changing  the  values  of 
intrinsic  permeability  and  displacement  pressure. 

Illangasekare  et  al.  (1995)  and  Kueper  and  Frind  (1991)  carried  out  DNAPL 
migration  studies  in  water-saturated  heterogeneous  porous  media.  These  cases  are  used 
for  comparison  to  the  stochastic  aggregation  model. 

5.2  Single  Fine/Coarse  Lens 

The  experiments  by  Illangasekare  et  al.  (1995)  used  a  flume  measuring  1.22  m 
v^dde  X  1.83  m  high  x  0.05  m  tiiick.  The  flume  walls  were  composed  ofV*  in.  thick 
plexiglas®  lined  with  Ya  in.  thick  tempered  glass  to  allow  visualization  of  the  experiment. 
A  steel  frame  supported  the  walls.  The  flume  was  packed  with  unconsolidated  sands  in 
two  configurations.  In  the  first  experiment,  a  20  cm  layer  of  #70  sand  was  placed  in 
otherwise  homogeneous  #30  sand.  In  the  second  experiment,  a  20  cm  layer  of  #16  sand 
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was  placed  in  the  homogeneous  #30  sand.  In  both  cases,  the  bottom  of  the  #30  or  #70 
layer  was  placed  40  cm  above  the  base  of  the  flume.  The  DNAPL  was  injected  100  cm 
above  the  base  of  the  flume  into  initially  saturated  sand.  Figure  5.1  from  Illangasekare  et 
al.  (1995)  shows  a  schematic  of  the  soil  configuration. 
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Figure  5.1.  The  experimental  configuration  used  by  Illangasekare  et  al.  (1995). 


1,1,1-trichloroethane  (TCA)  is  the  DNAPL  used  in  both  experiments.  TCA  has  a 
density  of  1.349  g  cm'^  and  a  dynamic  viscosity  of  1.2  cP  (Illangasekare  et  al.,  1995). 
Table  5.1  shows  the  displacement  pressure,  intrinsic  permeability,  transition  number  and 
sticking  probability  for  each  soil.  The  laboratory-scale  length  of  100  cm  and  an  injection 
rate  of  59  ml  min*'  are  used  in  the  calculation  of  the  transition  number. 


Table  5.1 

Sand  characteristics  and  model  parameters  for  the  experiments  of  Illangasekare  et  al.  (1995). 


Soil  type 

Intrinsic 

Permeability 

(cm^) 

DNAPL-water 

Displacement  Pressure 
(dyn/cm^) 

T 

Ps 

#16 

9.18E-6 

3588 

-1.58 

#30 

2.26E-6 

7791 

-3.09 

#70 

2.78E-7 

17932 

-11.1 

The  same  relationship  (Equation  [16])  is  used  to  relate  the  sticking  probability  to 


the  transition  number.  The  size  of  the  grid  corresponding  to  the  approximate  number  of 
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pores  in  the  flume  is  estimated  to  be  2844  pores  in  height  and  4220  pores  in  width.  The 
model  scaling  technique  discussed  in  Chapter  4  is  used  since  it  is  not  feasible  to  run 
simulations  of  this  size.  The  grid  size  for  the  simulation  is  the  same  as  the  grid  size  for 
the  model  validation,  113x113.  The  grid  must  be  square  to  use  the  graphics  program, 
PSPLOT,  to  visualize  the  output.  The  vertical  dimensions  of  the  flume  from  the  injection 
point  are  scaled  onto  a  1 13  x  1 13  grid  due  to  this  limitation  and  because  the  flume  was 
packed  synunetrically  about  the  vertical  centerline.  The  center  of  the  first  row  in  the  grid 
is  used  as  the  seed  point  corresponding  to  the  injection  point.  From  the  top  to  the  bottom 
of  the  flume:  the  first  40  cm  layer  of  #30  sand  corresponds  to  the  top  45  rows  of  the  grid, 
the  20  cm  layer  of  #70  sand  corresponds  to  the  next  23  rows  of  the  grid,  and  the  bottom 
40  cm  of  #30  corresponds  to  the  bottom  45  rows  in  the  grid. 

The  percentage  of  pore  space  occupied  by  the  DNAPL  is  estimated  to  compare 
simulations  with  the  physical  model  results  at  a  given  time.  The  volume  of  a  single  pore 
is  calculated  as  the  volume  of  a  sphere  with  the  radius  estimated  to  be  twice  the  radius  of 
a  pore  neck.  The  diameter  of  a  pore  neck  is  estimated  to  be  42%  of  dso.  the  average  grain 
diameter  (Ng  et  al.,  1978;  Pantazidou  and  Sitar,  1993;  Schroth  et  al.,  1995).  The  number 
of  pores  in  the  flume  is  approximated  by  taking  the  dimension  of  the  flume  in  each 
direction  and  dividing  it  by  grain  diameter  for  the  appropriate  sand  type.  The  volume  of  a 
single  pore  is  multiplied  by  the  total  number  of  pores  in  the  flume.  Multiplying  die 
number  of  pores  in  each  flume  dimension  and  summing  these  values  for  each  sand  type 
yields  an  approximate  pore  volume  of  28,600  cm^  for  the  experiment  with  the  #70  lens 
and  29,000  cm^  for  the  experiment  with  the  #16  lens.  The  percent  of  DNAPL-occupied 
pore  space  is  calculated  by  dividing  the  volume  of  DNAPL  by  the  volume  of  pore  space. 
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The  number  of  particles  necessary  for  a  simulation  is  determined  by  multiplying  the 
percentage  of  DNAPL  occupied  pore  space  by  the  number  of  grid  spaces  (1 13^).  Table 
5.2  summarizes  the  model  scaling  characteristics. 

Table  5.2 

Model  Scaling  Characteristics. 


Experiment  1 
#70  lens 

Experiment  2 
#16  lens 

Volume  ot  pore 
space  (cm^) 

28,600 

29,000 

%  of  TCA  occupied 
pore  space  at  ti 

3.11 

4.44 

Corresponding  number 

of  particles  at  ti 

397 

567 

%  of  TCA  occupied 

pore  space  at  t2 

9.32 

10.06 

Corresponding  number 

of  particles  at  t2 

1190 

1284 

Note:  ti  =  15  minutes  for  Experiment  1  and  30  minutes  for  Experiment  2 
t2  =  45  minutes  for  Experiment  1  and  68  minutes  for  Experiment  2 


Figure  5.2  illustrates  the  model  results  for  the  first  experiment.  Figure  5.2(a) 
compares  the  model  simulation  with  the  experiment  at  15  minutes.  Figures  5.2(b)  and 
5.2(c)  compare  two  different  model  simulations  with  the  experiment  15  hours  after  the 
spill.  In  Figure  5.2(b)  the  simulation  contains  the  volume  corresponding  to  total  volume 
injected.  It  is  more  appropriate  to  compare  this  simulation  with  a  photograph  taken  at  45 
minutes  (the  end  of  injection),  but  a  photograph  from  that  time  is  not  available.  The 
model  simulation  in  Figure  5.2(c)  shows  fluid  redistribution  and  spreading  at  the  textural 
interface  corresponding  to  the  photograph  taken  at  15  hours.  This  is  accomplished  by 
using  the  average  residual  saturation  of  0.3  measured  by  Illangasekare  et  al.  (1995)  to 
calculate  an  additional  number  of  particles  to  release.  The  area  occupied  by  TCA  in 
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The  model  simulations  qualitatively  match  the  experimental  pictures.  Both  the 
computer  simulations  and  the  physical  model  show  that  the  TCA  spreads  out  on  the  #70 
layer.  In  the  experiment,  TCA  is  not  able  to  penetrate  this  layer  because  the  displacement 
pressure  is  not  high  enough.  The  lowered  sticking  probability  accounts  for  this 
phenomenon  in  the  simulation.  Saturations  were  measured  by  Illangasekare  et  al.  (1995) 
using  gamma  attenuation  and  show  that  41 .6  hours  after  the  spill  the  TCA  saturations 
remain  nearly  zero  in  the  fine  layer.  Saturations  in  the  model  simulation  and  experiments 
are  not  easily  compared  as  the  gamma  attenuation  system  is  a  point  measurement  that  is 
not  meaningful  in  the  model  simulation. 

Figure  5.3  shows  the  results  of  the  model  simulation  of  the  second  experiment.  In 
this  experiment,  a  coarse  layer  of  #16  sand  is  sandwiched  between  layers  of  #30  sand. 
Figure  5.3(a)  is  the  DNAPL  configuration  after  the  first  30  minutes  of  injection.  Figures 
5.3(b)  and  5.3(c)  show  two  different  model  simulations  of  the  same  experiment 
photograph.  Figure  5.3(b)  shows  a  model  simulation  corresponding  to  the  time  at  which 
TCA  injection  stopped  (68  minutes).  Figure  5.3(c)  shows  a  model  simulation  where 
additional  particles  have  been  released  to  account  for  the  fluid  redistribution.  The 
experimental  photograph  in  Figures  5.3(b)  and  5.3(c)  was  taken  2  hours  after  the  spill. 

A  direct  comparison  between  the  model  simulations  and  experimental 
photographs  is  more  difficult  for  the  case  of  the  coarse  layer.  One  problem  is  that  the 
path  the  TCA  took  through  the  #16  layer  is  not  clearly  visible  in  the  photograph  in 
Figures  5.3(b)  and  5.3(c).  Illangasekare  et  al.  (1995)  indicate  that  fingers  were  observed 
in  the  #16  layer,  and  the  model  simulation  in  Figure  5.3(c)  shows  that  the  TCA 
configuration  within  the  #16  layer  is  more  fingered  compared  to  the  configuration  in  the 
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#30  layer.  The  model  simulation  in  Figure  5.3(c)  does  not  show  mounding  of  the  TCA 
on  the  lower  #30  layer. 


5.3(c) 

Figure  5.3.  Model  simulation  of  layered  system:  #30-#16-#30. 


In  the  first  experiment,  the  sticking  probability  corresponding  to  the  #70  sand  was 
0.005 1 .  This  is  small  enough  to  show  particles  collecting  before  eventually  sticking  in 
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this  region.  This  does  not  occur  in  the  model  simulation  of  the  coarse  layer  because  the 
sticking  probabilities  corresponding  to  the  #30  sand  and  the  #16  sand  are  relatively  large, 
0.0275  and  0.0669,  respectively.  An  evaluation  of  the  height  of  TCA  necessary  to  exceed 
the  entry  pressure  in  the  lower  #30  layer  is  missing.  This  model  limitation  is  discussed  in 
more  detail  in  Chapter  8.  The  longest  finger  will  continue  to  grow  preferentially  until  all 
of  the  particles  have  been  incorporated  into  the  aggregate  without  a  limitation  on  particle 
sticking  based  on  a  fluid  height. 

In  summary,  the  stochastic  aggregation  model  was  compared  to  two 
heterogeneous  laboratory  experiments  performed  by  Illangasekare  et  al.  (1995).  The 
simulations  were  performed  on  a  scaled  model  grid.  Both  the  grid  and  the  DNAPL 
volume  were  scaled.  This  was  done  by  approximating  the  percentage  of  pore  space 
occupied  by  the  DNAPL  and  allowing  the  number  of  particles  corresponding  to  this 
percentage  to  form  an  aggregate  onall3xll3  grid.  A  comparison  of  model 
simulations  and  laboratory  results  is  limited  to  visual  similarity  for  heterogeneous 
systems  (Kueper  and  Frind,  1991).  The  stochastic  aggregation  model  shows  good 
agreement  to  the  layered  experiments  of  Illangasekare  et  al.  (1995).  A  more  complex 
heterogeneous  system  is  considered  in  section  5.3. 
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5.3  Heterogeneous  Sand  Packing 

Model  simulations  are  compared  to  the  laboratory  experiment  performed  by 
Kueper  et  al.  (1989)  modeled  by  Kueper  and  Frind  (1991a  and  b).  In  their  experiment, 
tetrachloroethylene  (PERC)  was  allowed  to  infiltrate  into  a  flume  packed  with  four 
different  sands.  The  parallel-plate,  glass-lined  flume  was  60  cm  high  x  80  cm  wide  x  0.6 
cm  thick.  The  flume  was  packed  with  four  types  of  sands  in  the  configuration  shown  in 
Figure  5.4. 
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Figure  5.4.  The  experimental  configuration  for  Kueper  et  al.  (1989). 


PERC  was  injected  into  the  water-saturated  sand  at  the  top  of  the  flume.  The 
density  and  viscosity  of  PERC  is  1.63  g  cm'^  and  0.9  cP,  respectively.  The  top  and 
bottom  of  the  flume  were  sealed,  except  for  the  injection  area  shown  in  Figure  5.4.  The 
sides  of  the  flume  allowed  the  displaced  water  to  escape  through  a  fine  wire  screen. 
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The  DNAPL  infiltrated  into  the  water-saturated  flume  under  a  constant  head  of  4 
cm.  This  corresponds  to  a  transient  flow  rate  at  the  inflow  boundary.  A  volumetric  flow 
rate  is  estimated  for  each  of  six  PERC  distributions  corresponding  to  six  different  times 
shown  in  Figure  5.5.  The  flow  rate  is  estimated  by  enlarging  the  pictures  for  each 
distribution,  graphically  estimating  the  volume,  multiplying  the  volume  by  the  DNAPL 
saturation  given  by  Kueper  and  Frind  (1991a)  as  38%,  and  then  dividing  by  the  time  at 
which  the  distribution  is  observed.  The  flow  rates  do  not  vary  much  given  the  lack  of 
sophistication  of  the  method.  In  the  order  of  increasing  time  the  flow  rates  are  estimated 
to  be:  .Qa=0.586  cm^  s\  Qb=0.535  cm^  s’‘,  Qc=0.545  cm^  s'*,  Qd=0.591  cm^  s'*,  Qe=0.593 
cm^  s'*,  and  Qf=0.578  cm^  s'*. 

The  size  of  the  grid  corresponding  to  the  approximate  number  of  pores  in  the 
flume  is  estimated  to  be  424  vertically  and  593  horizontally.  This  is  based  on  the  size  of 
the  #16  sand.  The  model’s  grid  must  be  square  in  order  to  use  the  graphics  program  to 
visualize  the  output.  This  limitation  is  slightly  more  complex  in  this  case  because  the 
experimental  packing  of  the  flume  was  not  symmetric.  The  vertical  dimension  of  the  grid 
(424)  is  held  constant  and  the  horizontal  dimensions  are  scaled  by  71 .5%.  An  additional 
subroutine  is  implemented  into  the  stochastic  aggregation  model  containing  information 
on  the  grid  areas  corresponding  to  different  sticking  probabilities. 

The  number  of  particles  corresponding  to  the  volume  of  DNAPL  is  determined 
next.  The  given  dimensions  for  each  sand  type  are  multiplied  by  the  porosity  of  each 
sand  type  and  added  together  to  estimate  the  total  pore  volume.  The  total  pore  volume  is 
estimated  to  be  843  cm^  The  volume  of  DNAPL  is  divided  by  the  estimated  volume  of 
pore  space  to  give  a  percentage  of  pore  space  occupied  using  the  average  volumetric  flow 
rate  and  the  time.  The  percent  of  pore  space  occupied  is  multiplied  by  the  grid  size 
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(424^)  to  obtain  the  nvimber  of  particles  to  be  released  at  each  time  corresponding  to  the 
six  diagrams  shown  in  Figure  5.5.  Table  5.3  shows  the  estimated  volume  of  PERC, 
percent  of  occupied  pore  space,  and  corresponding  number  of  particles  for  each  time. 


Table  5.3 

Number  of  particles  for  model  simulation  of  Kueper  et  al.  (1989). 


Time  (s) 

Volume  of  PERC 
(cm') 

%  of  PERC 
occupied  pore  space 

Corresponding  number 
of  particles 

34 

19.42 

2.3 

4140 

126 

71.98 

8.5 

15350 

184 

105.11 

12.5 

22420 

220 

125.68 

14.9 

26800 

245 

139.96 

16.6 

29850 

313 

178.82 

21.2 

38130 

The  experiment  is  modeled  twice.  In  the  first  set  of  model  simulations,  the 
intrinsic  permeabilities  and  the  displacement  pressures  as  reported  by  Kueper  et  al. 

(1989)  for  each  sand  are  used.  Table  5.4  shows  the  soil  properties,  transition  number  and 


sticking  probability  for  each  sand  type  for  the  first  set  of  simulations.  The  laboratory- 


scale  length  of  60  cm  is  used  in  the  calculation  of  the  transition  number.  Figure  5.5 


shows  the  model  simulations  and  experimental  diagrams.  Photographs  of  the  experiment 


were  not  obtained,  but  would  be  preferable  for  comparison. 


Table  5.4 

Sand  characteristics  and  model  parameters  for  simulation  #1  of  Kueper  et  al.  (1989)  experiment. 


Soil  Type 

Intrinsic 

DNAPL-water 

T 

Ps 

Permeability 

Displacement  Pressure 

(cm') 

(dyn  cm'^) 

#16 

5E-6 

2266 

3.16 

#25 

2.1E-6 

2663 

6.41 

#50 

5.3E-7 

8116 

8.08 

0.00773 

#70 

8.2E-8 

19901 

21.05 

0.00217 
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Figure  5.5.  Model  simulations  compared  to  the  experimental  results  of  Kueper  et  al.  (1989) 


There  is  good  agreement  between  the  model  simulations  and  the  experiment.  The 
comparison  is  not  exact,  and  there  is  no  reason  to  believe  that  an  exact  comparison  is 
realistic.  The  model  simulations  shown  in  Figure  5.5  resulted  from  running  the  model 
once  for  the  given  set  of  inputs.  Likewise,  the  experimental  diagrams  show  the  result  of 
running  the  physical  experiment  once.  It  is  arguable  that  another  experimental  run  would 
not  produce  the  exact  same  results  even  with  attention  to  detail  in  the  packing  of  the 
flume  and  other  procedures  employed  in  the  experiment.  Results  refer  to  the  exact 
pathways  taken  by  the  DNAPL,  the  extent  to  which  the  DNAPL  travels  at  a  given  time, 
etc.  It  is  reasonable  that  the  same  experimental  setup  and  procedures  would  yield  similar 
results.  This  means  that  the  DNAPL’s  path  would  characteristically  favor  the  coarser 
sands  and  avoid  the  finer  sands.  This  is  also  true  of  the  model.  The  output  of  two  model 
simulations  run  with  the  same  number  of  particles,  corresponding  to  the  same  time  during 
the  experiment  are  shown  in  Figure  5.6.  They  show  characteristically  the  same  DNAPL 
migration  pattern.  They  are  not  exactly  the  same  but  show  that  as  the  DNAPL  moves,  it 
characteristically  avoids  the  finest  sand  and  preferentially  moves  through  the  coarser 
sands.  The  simulations  show  a  realistic  DNAPL  configuration. 


Figure  5.6.  Two  model  simulations  run  to  match  the  experiment  at  184  seconds  show 
characteristically  the  same  DNAPL  migration  path. 
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Kueper  and  Frind  (1991a)  modeled  the  experiment  of  Kueper  et  al.  (1989).  Their 
model  results  show  how  the  saturations  change  in  a  topographic  style  of  output.  This 
information  is  quantitatively  important,  but  cannot  visually  demonstrate  how  the  water 
and  DNAPL  occupy  the  pore  space. 

There  are  several  differences  between  the  model  and  experimental  results.  One 
difference  occurs  in  the  top  lens  marked  #3  in  Figure  5.4.  In  the  experiment,  the  PERC 
avoided  this  lens  and  moved  around  it.  The  model  does  not  capture  tiiis  behavior  for  two 
reasons.  First,  the  longest  finger  has  the  highest  probability  of  growth.  There  will  always 
be  particles  which  pass  through  the  lens  even  though  the  sticking  probability  for  the  #3 
sand  is  an  order  of  magnitude  lower  than  that  of  the  #1 .  As  time  progresses  more  and 
more  particles  find  their  way  around  the  #3  lens  (compare  Figure  5.5(c)  to  5.5(d)).  The 
second  reason  for  this  difference  is  that  the  model  does  not  check  the  entry  pressure  to 
determine  if  it  was  exceeded  for  the  #3  sand.  This  phenomenon  was  seen  in  the 
simulations  in  section  5.2.  Another  difference  between  the  experiment  and  the  model 
simulations  is  in  the  depth  that  the  PERC  travels.  The  model  simulation  shows  the  PERC 
traveling  farther  vertically  than  in  the  experiment.  The  model  simulations  do  not  show 
the  same  extent  of  lateral  spreading,  contributing  to  this  problem.  The  lack  of  lateral 
spreading  is  due,  in  part,  to  not  checking  the  entry  pressure.  If  the  model  evaluated  an 
entry  pressure  at  the  interface  between  the  #1  and  #3  sand,  particles  could  not  enter  the 
#3  sand  until  the  entry  pressure  is  exceeded.  This  would  lead  to  more  lateral  placement 
of  the  particles.  The  crude  estimation  of  the  grid  size  is  another  possible  explanation  for 
the  disparity  between  the  model  simulations  and  the  experiment.  The  grid  size  is 
estimated  by  assuming  that  there  is  the  same  number  of  pores  as  particles.  This 
assumption  relies  on  an  ideal  packing.  The  grid  size  is  calculated  based  on  only  one  sand 
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size,  the  #1  sand,  in  the  simulation  of  the  experiment  of  Kueper  et  al.  (1989).  The  #1  sand 
is  chosen  because  the  majority  of  the  DNAPL  traveled  through  the  #1  sand.  If  some 
combination  of  the  other  sand  sizes  is  used,  the  resulting  grid  size  would  be  larger 
because  the  grain  sizes  are  smaller.  Rvuming  the  simulations  on  a  larger  grid  would  result 
in  shallower  vertical  extent  of  the  DNAPL.  Estimation  of  the  grid  size  is  an  area  of 
fiiture  research. 

The  experiment  is  modeled  a  second  time  utilizing  the  tetrachloroethylene  source 
saturation  of  38%  used  by  Kueper  and  Frind  (1991a)  in  their  finite  difference  model.  The 
saturation  is  used  to  calculate  the  corresponding  capillary  pressure  in  the  #16  silica  sand 
using  the  Brooks-Corey  (1964)  capillary  presure-saturation  function  given  by: 


[17] 


where  Se  is  the  effective  saturation.  Pc  is  the  capillary  pressure,  Pd  is  the  displacement 
pressure,  and  X,  is  the  pore  size  distribution.  The  effective  saturation  can  be  calculated 
given  the  residual  saturation.  The  residual  saturation,  displacement  pressure,  and  pore 
size  distribution  are  given  by  Kueper  et  al.  (1989).  The  corresponding  capillary 
pressure,  calculated  to  be  5.03  cm  of  water,  is  used  in  the  transition  number  equation  for 
each  of  the  four  sands.  The  changes  in  the  transition  number  and  consequently  the 
sticking  probability,  are  due  solely  to  the  changes  in  grain  size  reflected  in  the  intrinsic 
permeability  values.  It  is  necessary  to  use  this  capillary  pressure  for  each  soil  because 
DNAPL  saturation  information  for  each  sand  is  not  known.  The  heterogeneities  for  this 
experiment  are  well  defined,  but  in  a  field  situation  the  heterogenities  could  be 
statistically  characterized  with  permeability  values.  Using  the  permeability  information. 
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the  saturation  at  the  injection  point,  and  information  on  the  statistical  distribution  of  the 
heterogenities,  the  transition  numbers  and  sticking  probabilities  could  be  calculated.  The 
values  of  the  transition  number  and  sticking  probability  for  each  sand  type  are  given  in 
Table  5.5.  The  model  simulations  and  experimental  diagrams  are  shown  in  Figure  5.7  for 
the  second  simulation. 


Table  5.5 

Sand  characteristics  and  model  parameters  for  simulation  #2  of  Kueper  et  al.  (1989)  experiment. 


Soil  Type 

Intrinsic 

DNAPL-water 

T 

Ps 

Permeability 

Displacement  Pressure 

(cm^) 

(dyn  cm'^) 

#16 

5E-6 

3027 

2.37 

#25 

2.1E-6 

3027 

5.65 

0.0124 

#50 

5.3E-7 

3027 

21.65 

0.00209 

#70 

8.2E-8 

3027 

138.4 

0.00018 

The  simulations  (Figure  5.7)  match  the  experimental  results  more  closely  than 
those  in  Figure  5.5.  The  sticking  probability  for  the  #3  sand  is  low  enough  to  discourage 
the  DNAPL  from  passing  through  the  lens.  As  a  result,  more  particles  travel  around  the 
lens,  filling  in  the  areas  where  the  #1  sand  is  located.  Increased  lateral  spreading  of  the 
DNAPL  reduced  the  depth  to  which  the  DNAPL  traveled.  The  DNAPL  in  the 
simulations  in  Figure  5.7  traveled  farther  than  in  the  physical  experiment.  This  is 
attributed  to  the  crude  grid  size  estimation.  It  is  important  to  remember  that  each 
simulation  is  separate  when  evaluating  the  simulation  pictures.  The  model  does  not  stop 
when  a  certain  number  of  particles  is  reached,  write  the  output  file,  and  then  continue. 
The  model  simulation  is  run  for  a  certain  number  of  particles  each  time.  The  simulations 
show  that  the  model  captures  the  characteristic  DNAPL  flow  path  in  the  same  way  each 
time. 
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In  summary,  the  stochastic  aggregation  model  was  compared  to  the  heterogeneous 
laboratory  experiment  performed  by  Kueper  et  al.  (1989).  The  grid  size  was  calculated 
based  on  the  particle  size  of  the  sand  used  in  the  experiment.  The  comparison  of  the 
model  simulations  and  the  experiment  are  limited  to  a  visual  assessment.  While  there  are 
some  differences  between  the  model  simulations  and  the  results  of  the  experiment,  the 
model  is  able  to  incorporate  a  complex  array  of  heterogenities  and  show  generally  the 
same  behavior.  Both  the  experiment  and  the  model  simulations  were  performed  only 
once.  It  uncertain  that  the  experiment  would  yield  the  same  output  even  if  performed  in 
exactly  the  same  way.  Likewise  the  model  simulations  show  only  one  simulation  for 
each  time.  Comparing  minute  details  of  the  simulations  with  the  experiment  does  not 
make  sense.  It  is  more  important  to  assess  whether  the  model  can  capture  the  general 
behavior  of  the  DNAPL  as  it  displaces  water  in  a  heterogeneous  porous  media. 
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CHAPTER  6 


MODEL  DETERMINATION  OF  BULK  RETENTION  CAPACITY 


Bulk  retention  capacity,  Rs,  is  defined  by  Johnson  and  Kueper  (1996)  as  the 
average  volume  of  NAPL  which  is  retained  per  unit  volume  of  the  aquifer.  The  bulk 
retention  capacity  differs  from  the  residual  saturation  of  a  NAPL  in  that  the  residual 
saturation  usually  refers  to  the  fi'action  of  pore  space  occupied  by  NAPL.  The  volume  of 
aquifer  includes  the  total  volume  of  soil,  gas,  and  liquid  through  which  the  DNAPL 
migrates.  According  to  the  AATDF  Technology  Practices  Manual  for  Surfactants  and 
Cosolvents:  Second  Edition  (1997),  the  volume  is  defined  as  the  “overall  volume  of 
medium  [and]  includes  both  those  lenses  and  laminations  in  which  residual  and  pooled 
NAPL  are  present  and  the  adjacent  lenses  and  laminations  void  of  NAPL”.  This  concept 
is  useful  in  providing  an  estimation  of  the  extent  of  the  contamination  in  the  field.  Most 
DNAPLs  have  extensive  vertical  migration  patterns  similar  to  the  fingers  shown  in  the 
laboratory  experiments  in  Chapter  3.  The  bulk  retention  capacity  for  a  given  volume  of 
soil  is  usually  quite  small  ranging  from  0.25  to  3  percent  (AATDF,  1997)  because  of 
these  fingers. 

The  dilemma  encountered  in  estimating  bulk  retention  capacities  lies  in  the 
volume  determination.  The  volume  as  defined  by  AATDF  (1997)  is  evaluated  on  a  site 
by  site  basis  and  cannot  be  known  with  certainty  without  excavating  the  DNAPL  and  soil 
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at  a  site.  It  is  useful  to  estimate  the  bulk  retention  capacity  using  the  stochastic 
aggregation  model.  Using  the  model  allows  a  visualization  of  the  occupied  pore  space, 
and  a  variety  of  “what  if’  scenarios  can  be  performed.  The  lab  setup  of  Illangasekare  et 
al.  (1995)  was  chosen  to  demonstrate  the  model’s  use  in  determining  bulk  retention 
capacities.  The  bulk  retention  capacities  can  be  compared  relatively.  The  effect  of 
different  soil  types  (#16,  #30,  and  #70  sand),  different  DNAPLs,  and  different  flow  rates 
are  investigated. 

6.1  Retention  Capacity  Determination 

The  model  simulations  in  Chapters  4  and  5  show  the  way  in  which  the  DNAPL 
displaces  the  water  and  moves  through  the  soil.  The  simulated  migration  patterns  are  the 
initial  characteristic  DNAPL  distribution.  After  the  source  of  DNAPL  is  stopped  there 
continues  to  be  DNAPL  movement  and  redistribution.  This  phenomenon  was  seen  by 
Illangasekare  et  al.  (1995)  in  the  experiments  simulated  in  Chapter  5. 

A  bulk  retention  capacity  for  a  specified  volume  is  evaluated  by  using  the  model 
to  determine  the  percentage  of  pore  space  affected  by  the  DNAPL  and  a  residual 
saturation.  Anderson  (1988)  reports  that  the  residual  saturation  for  PERC  in  saturated 
coarse  Ottawa  sand  ranges  from  0.15  to  0.25.  Mercer  and  Cohen  (1990)  report  the 
general  range  of  residual  saturations  for  DNAPLs  in  saturated  sand  are  0.02  to  0.15.  A 
residual  saturation  of  0.15  is  used  for  this  analysis.  Other  specified  parameters  include 
the  length  scale  of  100  cm  and  intrinsic  permeability  values  for  #16,  #30,  and  #70  sand  as 
measured  by  Illangasekare  et  al.  (1995).  The  bulk  retention  capacity  is  reported  as  L  m'^. 
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Table  6.1  contains  the  bulk  retention  capacity  calculated  for  5  common  DNAPLs 
in  three  different  homogeneous  sands  at  an  injection  rate  of  59  ml  min’*.  Tables  6.2  and 
6.3  contain  the  bulk  retention  capacities  for  the  flow  rates  of  30  ml  min  *  and  15  ml  min  \ 
respectively.  PERC,  CTET,  and  TCE  are  evaluated  because  they  are  common 
chlorinated  solvents.  Creosote  is  also  a  common  DNAPL  and  was  chosen  for  its  high 
viscosity.  Bromoform  was  chosen  for  this  example  because  it  is  more  viscous  than  water, 
but  not  as  viscous  as  creosote.  The  bulk  retention  capacities  range  from  0.62  to  2.86 
percent  by  bulk  volume.  This  is  close  to  the  same  range  reported  by  the  AATDF  (1997) 
of  0.25  to  3  percent  by  bulk  voliune. 


Table  6.1 

Model  calculated  retention  capacity  values  for  Q=59  ml  min'*. 


Sand 

DNAPL 

Density 

(g  cm'O 

Viscosity 

(cP) 

Area  Ratio 

%DNAPL/ 
Bulk  Volume 

Retention  Capacity 
(Lm-^) 

#16 

PERC 

1.63 

0.932 

0.3156 

1.467 

14.67 

CTET 

1.58 

0.965 

0.2852 

1.326 

13.26 

TCE 

1.47 

0.560 

0.4760 

2.213 

22.13 

Creosote 

1.03 

54.0 

0.6157 

2.863 

28.63 

Bromoform 

2.9 

2.15 

0.5349 

2.488 

24.88 

#30 

PERC 

1.63 

0.932 

0.4113 

1.913 

19.13 

CTET 

1.58 

0.965 

0.3813 

1.773 

17.73 

TCE 

1.47 

0.560 

0.5383 

2.503 

25.03 

Creosote 

1.03 

54.0 

0.6162 

2.865 

28.65 

Bromoform 

2.9 

2.15 

0.5751 

2.674 

26.74 

#70 

PERC 

1.63 

0.932 

0.5507 

2.561 

25.61 

CTET 

1.58 

0.965 

0.5365 

2.494 

24.94 

TCE 

1.47 

0.560 

0.5966 

2.774 

27.74 

Creosote 

1.03 

54.0 

0.6165 

2.867 

28.67 

Bromoform 

2.9 

2.15 

0.6066 

2.821 

28.21 
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Table  6.2  ^ 

Model  calculated  retention  capacity  values  for  Q=30  ml  min**. 


Sand 

DNAPL 

Density 
(g  cm'O 

Viscosity 

(cP) 

Area  Ratio 

%DNAPL/ 
Bulk  Volume 

Retention  Capacity 
(Lm-') 

#16 

PERC 

1.63 

0.932 

0.2186 

1.016 

10.16 

CTET 

1.58 

0.965 

0.1948 

0.906 

9.06 

TCE 

1.47 

0.560 

0.3799 

1.766 

17.66 

Creosote 

1.03 

54.0 

0.6146 

2.858 

28.58 

Bromoform 

2.9 

2.15 

0.4601 

2.139 

21.39 

#30 

PERC 

1.63 

0.932 

0.3039 

1.413 

14.13 

CTET 

1.58 

0.965 

0.2735 

1.272 

12.72 

TCE 

1.47 

0.560 

0.4690 

2.181 

21.81 

Creosote 

1.03 

54.0 

0.6156 

2.863 

28.63 

Bromoform 

2.9 

2.15 

0.5318 

2.473 

24.73 

#70 

PERC 

1.63 

0.932 

0.4894 

2.275 

22.75 

CTET 

1.58 

0.965 

0.4658 

2.166 

21.66 

TCE 

1.47 

0.560 

0.5742 

2.669 

26.69 

Creosote 

1.03 

54.0 

0.6164 

2.866 

28.66 

Bromoform 

2.9 

2.15 

0.5950 

2.767 

27.67 

Table  6.3 

Model  calculated  retention  capacity  values  for  Q=15  ml  min*‘. 

Sand 

DNAPL 

Density 
(g  cm'O 

Viscosity 

(cP) 

Area  Ratio 

%DNAPL/ 
Bulk  Volume 

Retention  Capacity 
(Lm-^) 

#16 

PERC 

1.63 

0.932 

0.2186 

0.691 

6.91 

CTET 

1.58 

0.965 

0.1948 

0.621 

6.21 

TCE 

1.47 

0.560 

0.3799 

1.260 

12.60 

Creosote 

1.03 

54.0 

0.6121 

2.846 

28.46 

Bromoform 

2.9 

2.15 

0.4601 

1.614 

16.14 

#30 

PERC 

1.63 

0.932 

0.3039 

0.950 

9.50 

CTET 

1.58 

0.965 

0.2735 

0.844 

8.44 

TCE 

1.47 

0.560 

0.4690 

1.708 

17.08 

Creosote 

1.03 

54.0 

0.6144 

2.857 

28.57 

Bromoform 

2.9 

2.15 

0.5318 

2.113 

21.13 

#70 

PERC 

1.63 

0.932 

0.4894 

1.831 

18.31 

CTET 

1.58 

0.965 

0.4658 

1.688 

16.88 

TCE 

1.47 

0.560 

0.5742 

2.461 

24.61 

Creosote 

1.03 

54.0 

0.6161 

2.865 

28.65 

Bromoform 

2.9 

2.15 

0.5950 

2.651 

26.51 
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6.2  Changes  in  Bulk  Retention  Capacity  with  Varying  Parameters 

The  bulk  retention  capacities  calculated  by  the  model  are  investigated  with  the 
variations  in  flow  rate,  intrinsic  permeability,  and  viscosity.  Figure  6.1  shows  the  bulk 
retention  capacity  increasing  as  the  flow  rate  increases  for  #30  sand.  As  the  flow  rate 
increases,  the  DNAPL  configuration  progresses  from  single  fingers  to  multiple  fingers 
and  the  DNAPL  occupies  more  of  the  area.  The  model  accounts  for  the  increased 
amoimt  of  area  occupied  by  the  DNAPL  without  generalizing  the  area  due  to  the  increase 
in  lateral  spreading  of  the  DNAPL.  The  retention  capacity  for  creosote  does  not  show  as 
dramatic  an  increase  as  the  other  DNAPLs,  due  to  its  high  viscosity.  The  data  for  each 
DNAPL  in  Figure  6.1,  including  creosote,  can  be  modeled  using  varying  coefficients  of  a 
simple  logarithmic  equation,  Rs  =  a  +  b  ln(Q),  where  a  and  b  are  empirical  coefficients 
determined  for  each  fluid. 
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Figure  6.1.  The  bulk  retention  capacity  increases  as  the  flow  rate 
increases  for  6  common  DNAPLs. 
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Figure  6.2  shows  the  change  in  bulk  retention  capacity  with  the  change  in  intrinsic 
permeability  for  a  flow  rate  of  30  ml  s'*.  The  other  flow  rates  show  similar  trends.  The 
bulk  retention  capacity  decreases  as  the  intrinsic  permeability  increases.  As  fluids  invade 
finer  sands,  the  capillary  forces  have  increasing  influence  over  fluid  movement,  and  the 
fluid  spreads  laterally  leading  to  a  greater  amount  of  entrapment  (Morrow  and  Songkran, 
1981).  The  amoimt  of  entrapped  fluid  (residual  saturation)  remains  the  same  (15%)  for 
this  example.  Figure  6.2  shows  the  influence  of  the  increased  spreading  due  to  the 
changing  intrinsic  permeability.  The  data  points  for  creosote  appear  to  form  a  straight 
line,  but  in  fact  follow  the  same  trend  as  the  other  DNAPLs.  The  curves  in  Figure  6.2  for 
each  DNAPL  can  be  modeled  with  a  simple  power  fit,  Rs  =  a  (k)** ,  with  empirical 
coefficients,  a  and  b  changing  for  each  DNAPL. 
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Figure  6.2.  The  bulk  retention  capacity  decreases  as  the  intrinsic 
permeability  increases  for  6  common  DNAPLs. 


67 


Figures  6.3  and  6.4  show  the  relationship  between  the  bulk  retention  capacity  and 
viscosity  difference  for  each  type  of  sand.  The  viscosity  difference  between  the  DNAPL 
and  water  increases  as  the  bulk  retention  capacity  increases  until  the  bulk  retention 
capacity  approaches  29  L  m'^.  In  Figure  6.4,  the  data  points  corresponding  to  creosote 
are  excluded  to  show  the  increasing  part  of  the  curve.  Figures  6.3  and  6.4  show  that 
fluids  with  high  viscosity  differences  relative  to  water,  such  as  creosote,  have  higher 
retention  capacities.  The  influence  of  different  intrinsic  permeabilities  is  also  seen  in 
these  figures.  The  values  corresponding  to  the  finest  sand  (#70)  have  (he  highest 
retention  capacities  for  the  same  DNAPL.  The  data  can  be  modeled  by  a  modified 
exponential  function  of  the  form,  Rs=a(b-e‘®^),  where  a,  b,  and  c  are  empirical  coefficients 
determined  for  each  sand  type. 
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Figure  6.3.  The  bulk  retention  capacity  approaches  29  L  m""*  with 
increasing  viscosity  difference  relative  to  water. 
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Figure  6.4.  The  portion  of  Figure  6.3  with  viscosity  differences  ranging  between  0  and  2. 


Changes  in  bulk  retention  capacity  with  changes  in  flow  rate,  intrinsic 
permeability,  and  viscosity  are  easily  assessed  by  using  the  model.  As  the  flow  rate  is 
increased  the  bulk  retention  capacity  increases  logarithmically.  This  is  attributed  to  a 
higher  capillary  pressme  gradient  yielding  more  pores  accessible  for  displacement.  The 
bulk  retention  capacity  increases  according  to  a  modified  exponential  function  as  the 
viscosity  difference  increases.  A  limit  of  28.7  L  m'^  is  approached  for  the  bulk  retention 
capacity  in  both  cases.  This  is  an  artifact  of  the  bulk  volume  chosen.  If  the  chosen  bulk 
volume  is  greater,  the  corresponding  model  grid  will  be  bigger  and  die  retention  capacity 
reaches  a  smaller  limit.  The  DNAPL  remains  in  2.87%  of  the  bulk  volume  in  the  limit  of 
very  high  flow  rate  or  very  large  viscosity  difference.  This  corresponds  to  the  DNAPL 
traveling  through  61%  of  the  pore  space  with  a  residual  saturation  of  15%.  Figure  6.5 
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shows  a  model  simulation  where  the  DNAPL  injected  from  a  point  source,  traveled 
through  approximately  61%  of  the  pore  space  in  the  given  volume. 


Figure  6.5.  A  DNAPL  occupying  61%  of  the  pore  space. 

The  model  is  also  used  to  investigate  the  effect  of  changing  intrinsic  permeability 

on  the  bulk  retention  capacity.  In  this  case  die  bulk  retention  capacity  decreased  with 
increasing  permeability  reaching  a  limit  dependent  on  the  DNAPL’s  viscosity  and  density 
properties. 

In  summary,  changes  in  bulk  retention  capacity  with  changing  fluid  and  soil 
properties  are  useful  for  field  scale  determinations.  The  stochastic  aggregation  model  can 
be  used  for  this  purpose.  The  bulk  retention  capacity  can  be  approximated  and  a 
visualization  of  the  way  in  which  the  DNAPL  initially  occupied  the  soil  is  useful  in 
characterizing  the  source  zone.  Varying  parameters  to  determine  a  range  of  values 
corresponding  to  differing  scenarios  due  to  die  imcertainty  in  parameter  values  is  easily 
accomplished. 
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CHAPTER  7 


MODEL  EVALUATION 

The  benefits  and  limitations  of  the  stochastic  aggregation  modeling  approach  are 
explored  in  this  chapter. 

7.1  Benefits  of  the  Stochastic  Aggregation  Model 

The  stochastic  aggregation  model  can  visually  simulate  the  characteristic  unstable 
flow  patterns  seen  in  DNAPL-water  displacement  in  porous  media.  The  objectives 
defined  during  the  formulation  of  the  model  included  linking  the  DLA  aggregate  growth 
to  the  front  instability,  exploration  of  the  model  to  show  trends  in  DNAPL  migration  with 
changing  fluid  properties  and  buoyancy,  incorporation  of  porous  media  heterogenities, 
determination  of  bulk  retention  capacity  values,  and  comparison  with  lab  experiments. 
These  objectives  were  satisfied.  Further  exploration  of  the  model  will  satisfy  larger  and 
arguably  more  relevant  field  scale  objectives. 

The  stochastic  aggregation  model  allows  the  user  to  visualize  the  migration 
pattern  of  a  DNAPL  traveling  through  water-saturated  porous  media  at  the  pore  scale. 

The  migration  pattern  is  based  on  the  porous  media  and  DNAPL  properties  including 
intrinsic  permeability,  displacement  pressure,  inflow  rate,  DNAPL  density  and  DNAPL 
viscosity.  The  model  does  not  generate  a  unique  DNAPL  configuration  for  a  given  set  of 
inputs  because  random  numbers  are  used  to  determine  the  placement  of  particles  on  the 
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grid.  The  model  simulations  show  characteristic  DNAPL  configurations.  The  important 
information  fi-om  the  model  realizations  is  the  characteristic  DNAPL  migration  pattern 
(fingered  or  flat  firont)  and  the  relative  amount  of  pore  space  occupied  as  the  DNAPL 
migrates  through  the  soil.  The  relative  amount  of  pore  space  occupied  by  the  DNAPL  is 
correlated  to  the  model  inputs  for  a  given  grid  size  and  remains  constant  for  fiiose  inputs. 
More  realistic  source  zone  configurations  are  the  ultimate  goal  of  the  model  simulations. 

The  stochastic  aggregation  algorithm  is  simple.  Complex  numerical  modeling 
techniques  are  not  needed.  Model  simulations  run  in  less  than  15  minutes  for  sticking 
probabilities  as  small  as  1  x  10'^  at  a  grid  size  of  113x113.  The  time  necessary  to  run 
simulations  increases  dramatically  vrith  increasing  grid  size.  It  was  shown  that  over  an 
order  magnitude  increase  in  grid  size,  the  characteristics  of  the  aggregate  on  a  given  scale 
remain  the  same,  and  allow  a  smaller  grid  to  be  used. 

Heterogenities  are  easily  incorporated  into  the  model.  Different  areas  of  the  grid 
are  assigned  different  sticking  probabilities  corresponding  to  changes  in  the  soil 
parameters.  Comparisons  between  model  simulations  and  the  experiments  of 
Illangasekare  et  al.  (1995)  and  Kueper  and  Frind  (1991)  show  good  agreement.  This  is 
the  first  study  of  a  stochastic  aggregation  algorithm  based  on  DLA  that  is  compared  to 
physical  experiments. 

The  model  is  used  to  calculate  bulk  retention  capacities.  Results  of  simulations  in 
this  study  are  within  the  range  of  retention  capacities  seen  in  field  experiments  (Poulsen 
and  Kueper,  1992)  and  cited  by  others  (AATDF,  1997).  The  value  of  the  model  is  that 
in  addition  to  a  relevant  bulk  retention  capacity,  a  graphic  of  the  DNAPL  migration  is 
also  available.  Using  the  model  to  determine  bulk  retention  capacities  allows  for  “what 
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if’  scenarios  to  be  run  in  field  situations  where  parameters  may  not  be  known  with 
certainty.  Other  methods  of  calculating  and/or  measuring  bulk  retention  capacities  in  tiie 
field  include  the  use  of  partitioning  tracers.  Partitioning  tracers  yield  a  field-scale  bulk 
average  NAPL  saturation.  It  is  possible  to  use  this  field-measured  value  to  back-calculate 
a  sticking  probability.  The  sticking  probability  may  then  be  used  as  model  input  to  create 
a  visualization  of  the  characteristic  DNAPL-water  displacement. 

7.2  Limitations  of  the  Stochastic  Aggregation  Model 

Many  limitations  of  the  stochastic  aggregation  model  have  been  discovered 
during  the  course  of  this  research.  This  is  expected  when  using  a  relatively  new 
technique.  Some  limitations  are  interesting  in  that  they  are  unique  to  this  modeling 
technique. 

The  lack  of  theory  between  DLA-type  models  and  fluid-fluid  displacement 
remains  the  largest  detriment  to  the  modeling  technique.  The  formulation  of  the  model  is 
based  on  a  calibration  of  the  numerical  description  of  the  stability  of  the  front  (T),  as  it  is 
influenced  by  the  fluid  and  porous  media  properties,  to  the  sticking  probability  in  the 
model.  In  this  first  attempt,  the  author  calibrated  the  model  using  a  set  of  laboratory 
experiments  instead  of  taking  on  the  problem  of  establishing  the  theoretical  link.  The 
general  lack  of  other  DNAPL  experiments  and  their  varying  laboratory  procedures  made 
broad  model  calibration  difficult.  The  lack  of  theory  pertains  to  the  relationship  between 
the  sticking  probability  and  the  physical  parameters  governing  the  movement  of  the 
DNAPL  as  it  displaces  water.  This  link  was  investigated  in  the  past  in  the  absence  of 
gravity  (Fernandez  et  al.,  1991;  Kiriakidis  et  al.,  1991)  but  has  not  before  been  linked  to 
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the  Saffman  and  Taylor  (1958)  and  Chuoke  et  al.  (1959)  equation  describing  the 
instability  of  the  front  Until  the  mathematics  describing  the  aggregate  growth  as  a 
function  of  a  changing  sticking  probability  are  determined,  the  theoretical  link  cannot  be 
established. 

Time  is  not  incorporated  into  this  stochastic  aggregation  model.  The  model  is 
only  applicable  to  steady  flow  conditions  at  the  inflow  boundary.  The  model  simulates 
DNAPL-water  displacement  at  varying  times  if  the  volume  of  DNAPL  present  at  the  time 
in  question  is  known.  If  the  volumetric  flow  rate  does  not  change  substantially  there  is 
good  agreement  between  the  model  simulation  and  the  experiment.  Other  researchers 
(Flxuy  and  Fluhler,  1995)  have  devised  time  aspects  into  modified  DLA  algorithms.  This 
is  an  area  for  future  expansion  of  the  stochastic  aggregation  model. 

The  model  does  not  determine  the  capillary  pressure  during  a  model  simulation 
when  varying  soil  properties  are  incorporated  into  the  grid.  The  displacement  pressure  is 
used  in  the  determination  of  the  transition  number,  and  it  provides  general  information 
about  the  porous  media.  In  this  capacity  it  does  not  ensure  that  the  aggregate  generated 
by  the  model  will  grow  in  such  a  way  that  areas  of  the  porous  media  having  high  entry 
pressures  will  not  be  invaded.  An  entry  pressure  incorporated  into  the  model  would 
allow  the  model  simulations  to  be  more  realistic.  The  aggregate  would  tend  to  grow 
preferentially  along  (lateral  to)  a  layer  of  lower  permeability  by  applying  the  entry 
pressure  as  a  constraint  on  the  addition  of  particles  to  the  aggregate.  The  same  volume  of 
DNAPL  would  not  travel  as  deep  with  more  lateral  spreading. 

There  are  cases  of  DNAPL-water  displacement  that  are  beyond  the  limits  of  the 
model  as  determined  by  tiie  model  calibration  and  cannot  be  simulated.  Combinations  of 
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parameters  resulting  in  transition  number  values  lower  than  0.205141  result  in  a  sticking 
probability  greater  than  one.  One  is  the  largest  sticking  probability  that  can  be  used  in 
the  model.  Each  sticking  probability  corresponds  to  an  area  ratio  given  by  Equation  [15]. 
The  cases  where  the  sticking  probability  is  greater  than  one  result  in  configurations  where 
the  area  ratio  is  less  than  the  smallest  sized  finger  the  model  can  generate.  The  model 
validation  cases  of  CTET  and  PERC  at  15®  are  included  to  show  that  using  a  sticking 
probability  of  one  works  reasonably  well  in  these  cases. 

Simulations  on  very  large  grids  (larger  than  1000  x  1000)  are  not  feasible  due  to 
large  computational  time  requirements.  A  simulation  run  with  a  sticking  probability  of 
0.8  took  22  hours  of  real  time  and  730  minutes  of  CPU  time  (IBM  RS6000)  on  a  1000  x 
1000  grid.  There  are  some  interesting  questions  that  could  be  investigated  if  it  were 
easier  to  run  larger  problems.  One  such  question  was  encountered  during  the  evaluation 
of  the  model  scaling.  It  is  not  clear  whedier  or  not  the  area  ratio  determined  for  a  set  grid 
size  and  a  given  sticking  probability  increases  with  grid  sizes  larger  than  1000  x  1000. 
This  is  an  important  consideration  if  the  model  is  to  be  used  on  field  scale  cases  where 
there  could  conceivably  be  grid  sizes  of  10^  x  10^  or  larger. 

The  current  model  performs  two-dimensional  simulations.  It  is  straight  forward 
to  design  a  similar  algorithm  in  three  dimensions.  A  three-dimensional  algorithm 
requires  a  theoretical  link  between  the  sticking  probability  and  the  physical  parameters. 
This  is  due  to  a  lack  of  three-dimensional  experiments  available  for  calibrating  the  model. 

The  graphics  program  used  to  generate  a  post  script  file  for  the  model’s  output 
limits  the  visualization  of  the  output.  This  is  not  a  limitation  with  the  model  itself,  but  of 
the  program  entitled  PSPLOT.  PSPLOT  was  obtained  with  written  permission  from 
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Kevin  Kohler,  Coordinator  of  Computing  Services  at  Nova  Southeastern  University 
Oceanographic  Center,  Dania,  FL.  The  code  is  written  in  Fortran  90  and  is  a  powerful 
and  extensive  program.  PSPLOT  shows  the  model’s  output  on  a  square  grid  of  specified 
size.  The  size  of  the  picture  ouqiut  by  PSPLOT  remains  the  same  regardless  of  the  grid 
size.  This  is  a  drawback  because  more  and  more  lines  are  drawn  in  the  same  amount  of 
space  as  the  grid  size  increases.  The  pictures  become  darker  and  darker.  The  resulting 
graphic  is  completely  black  at  a  grid  size  of  800  x  800. 

This  modeling  technique  shows  promise.  This  is  evidenced  by  the  close  matching 
of  the  model  simulations  with  the  laboratory  experiments.  This  model,  however,  is  new 
and  the  limitations  outlined  above  need  to  be  addressed. 
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CHAPTERS 


SUMMARY,  CONCLUSIONS,  AND  RECOMMENDATIONS 

The  objectives  of  this  research  were  to  model  DNAPL  migration  in  saturated 
porous  media  at  the  pore  scale,  evaluate  the  modeling  technique  through  comparisons 
with  DNAPL-water  displacement  studies  in  homogeneous  and  heterogeneous  porous 
media,  and  to  investigate  the  application  of  the  model  to  the  evaluation  of  bulk  retention 
capacities.  The  interface  between  the  fluids  is  often  unstable  due  to  DNAPL  properties. 
Analytical  solutions  stemming  from  the  continuity  equation  do  not  explicitly  consider  the 
stability  of  the  interface  where  displacement  occurs  (at  the  pore  scale).  Consequently, 
conventional  models  have  not  achieved  DNAPL  migration  and  source  zone 
characterization  at  the  pore  scale.  A  different  modeling  technique  similar  to  diffusion 
limited  aggregation  (DLA),  developed  by  Witten  and  Sander  (1983),  was  investigated  for 
this  purpose.  A  miique  relationship  between  the  model  input  parameter,  the  sticking 
probability,  and  the  essential  properties  governing  the  DNAPL-water  front  displacement 
was  developed.  The  stochastic  aggregation  model  is  the  first  DLA-based  model 
successfully  used  to  simulate  a  variety  of  laboratory  experiments,  both  in  homogeneous 
and  heterogeneous  porous  media. 

An  evaluation  of  the  model  concludes  that  the  model  is  well  formulated  for 
simulating  two-dimensional  DNAPL-water  displacement  for  a  variety  of  DNAPLs  under 
the  conditions  of  different  inflow  rates,  different  angles  of  inclination  from  the  horizontal. 
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and  varying  heterogeneous  media  properties.  The  comparison  of  model  simulations  to 
laboratory  experiments  remains  limited  to  visual  inspection.  This  is  the  case  for  other 
models  currently  attempting  to  model  DNAPL-water  displacements. 

The  model  was  used  to  evaluate  the  bulk  retention  capacity  for  a  variety  of 
common  DNAPLs.  A  bulk  retention  capacity  can  be  determined  for  a  specified  bulk 
volume  of  soil,  and  the  characteristic  front  displacement  can  be  visualized  using  the 
stochastic  aggregation  model.  The  bulk  retention  capacity  increases  as  the  flow  rate  of 
the  DNAPL  increases,  as  the  viscosity  difference  with  water  increases,  and  as  the 
intrinsic  permeability  of  the  media  decreases. 

The  use  of  this  modeling  technique  in  evaluating  the  problem  of  DNAPL  source 
zone  characterization  is  new.  It  is  not  yet  developed  to  a  stage  that  is  ready  for  broad 
field-scale  usage.  There  are  many  questions  remaining  to  be  answered  before  a  broad 
application  and  extension  of  the  model  can  occur.  The  following  are  recommendations 
for  future  research. 

•  The  theoretical  relationship  between  the  capillary  pressure  gradient  and  the 
model’s  sticking  probability  should  be  established.  This  caimot  be  done  until  the 
mathematics  behind  incorporating  a  sticking  probability  into  the  model  are  known.  If  the 
theoretical  link  can  be  established,  the  model  can  be  used  to  simulate  three-dimensional 
displacement  fronts. 

•  Scaling  of  the  grid  size  should  be  investigated  further  by  running  very  large 
simulations.  The  pmpose  of  scaling  the  grid  is  to  reduce  the  time  necessary  for 
simulations.  It  is  reasonable  to  believe  that  the  advancement  of  technology  may  take  care 
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of  this  problem.  There  may  no  longer  be  a  need  for  grid  scaling  if  the  time  it  takes  to  run 
large-sized  grid  simulations  is  improved. 

•  Grids  of  stochastically  varying  permeabilities  should  be  generated  and  used 
with  the  model  to  produce  simulations  that  are  more  representative  of  field  conditions. 

•  An  evaluation  of  the  entry  pressure  should  be  incorporated  into  the  model. 
Incorporating  an  entry  pressure  determination  into  the  model  gives  the  model  another 
piece  of  information  on  which  to  base  decisions  for  the  firont  movement.  This  entry 
pressure  evaluation  will  lead  to  more  accurate  model  simulations. 
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